Fully Discrete Finite Element Approximation for the Stabilized Gauge-Uzawa Method to Solve the Boussinesq Equations

Abstract

The stabilized Gauge-Uzawa method (SGUM), which is a 2nd-order projection type algorithm used to solve Navier-Stokes equations, has been newly constructed in the work of Pyo, 2013. In this paper, we apply the SGUM to the evolution Boussinesq equations, which model the thermal driven motion of incompressible fluids. We prove that SGUM is unconditionally stable, and we perform error estimations on the fully discrete finite element space via variational approach for the velocity, pressure, and temperature, the three physical unknowns. We conclude with numerical tests to check accuracy and physically relevant numerical simulations, the Bénard convection problem and the thermal driven cavity flow.

1. Introduction

The stabilized Gauge-Uzawa method (SGUM) is a 2nd-order projection type method to solve the evolution Navier-Stokes equations. In this paper, we extend SGUM to the evolution Boussinesq equations: given a bounded polygon Ω in ℝ^𝑑 with 𝑑=2 or 3, 𝐮_𝑡+(𝐮⋅∇)𝐮+∇𝑝−𝜇Δ𝐮+𝜅𝜇²𝜃𝐠𝜃=𝐟,inΩ,∇⋅𝐮=0,inΩ,_𝑡+𝐮⋅∇𝜃−𝜆𝜇Δ𝜃=𝑏,inΩ, (1) with initial conditions 𝐮(𝐱,0)=𝐮⁰,  𝜃(𝐱,0)=𝜃⁰ in Ω, vanishing Dirichlet boundary conditions 𝐮=0,  𝜃=0 on 𝜕Ω, and pressure mean-value ∫_Ω𝑝=0. The forcing functions 𝐟 and 𝑏 are given, and 𝐠 is the vector of gravitational acceleration. The nondimensional numbers 𝜇=Re⁻¹ and 𝜆=Pr⁻¹ are reciprocal of the Reynolds and Prandtl numbers, respectively, whereas 𝜅 is the Grashof number. The Boussinesq system (1) describes fluid motion due to density differences which are in turn induced by temperature gradients: hot and thus less dense fluid tends to rise against gravity and cooler fluid falls in its place. The simplest governing equations are thus the Navier-Stokes equations for motion of an incompressible fluid, with forcing 𝜅𝜇²𝐠𝜃 due to buoyancy and the heat equation for diffusion and transport of heat. Density differences are thus ignored altogether except for buoyancy.

The projection type methods are representative solvers for the incompressible flows, and the Gauge-Uzawa method is a typical projection method. The Gauge-Uzawa method was constructed in [1] to solve Navier-Stokes equations and extended to more complicated problems, the Boussinesq equations in [2] and the nonconstant density Navier-Stokes equations in [3]. However, most of studies for the Gauge-Uzawa method have been limited only for the first-order accuracy backward Euler time marching algorithm. The second-order Gauge-Uzawa method using BDF2 scheme was introduced in [4] and proved superiority for accuracy on the normal mode space, but we could not get any theoretical proof via energy estimate even stability and we suffer from weak stability performance on the numerical test. Recently, we construct SGUM in [5] which is unconditionally stable for semidiscrete level to solve the Navier-Stokes equations. The goal of this paper is to extend SGUM to the Boussinesq equations (1), which model the motion of an incompressible viscous fluid due to thermal effects [6, 7]. We will estimate errors and stability on the fully discrete finite element space. The main difficulties in the fully discrete estimation arise from losing the cancellation law due to the failing of the divergence free condition of the discrete velocity function. The strategy of projection type methods computes first an artificial velocity and then decomposes it to divergence free velocity and curl free functions. However, the divergence free condition cannot be preserved in discrete finite element space, and so the cancellation law (12) can not be satisfied any more. In order to solve this difficulty, we impose the discontinuous velocity on across interelement boundaries to make fulfill discrete divergence free velocity (12) automatically. We will discuss this issue at Remark 2 below. This discontinuity makes it difficult to treat nonlinear term and to apply the integration by parts, because the discontinuous solution is not included in 𝐻¹(Ω). So we need to hire technical skills in proof of this paper.

One more remarkable discovery is in the second numerical test at the last section which is the Bénard convection problem with the same setting in [2]. In this performance, we newly find out that the number of circulations depends on the time step size 𝜏. We obtain similar simulation within [2] for a relatively larger 𝜏, but the behavior is changed for the small 𝜏. So we conclude that the numerical result of the Bénard convection problem in [2] is not an eventual solution and a more smaller time marching step is required to get the desired simulation.

We will denote 𝜏 as the time marching size. Also we will use 𝛿 as the difference of 2 consecutive functions for example, for any sequence function 𝑧^𝑛+1, 𝛿𝑧^𝑛+1=𝑧^𝑛+1−𝑧^𝑛,𝛿𝛿𝑧^𝑛+1=𝛿𝛿𝑧^𝑛+1=𝑧^𝑛+1−2𝑧^𝑛+𝑧^𝑛−1,…. (2)

In order to introduce the finite element discretization, we need further notations. Let 𝐻^𝑠(Ω) be the Sobolev space with 𝑠 derivatives in 𝐿²(Ω), set 𝐋²(Ω)=(𝐿²(Ω))^𝑑 and 𝐇^𝑠(Ω)=(𝐻^𝑠(Ω))^𝑑, where 𝑑=2 or 3, and denote by 𝐿²₀(Ω) the subspace of 𝐿²(Ω) of functions with vanishing mean value. We indicate with ‖⋅‖_𝑠 the norm in 𝐻^𝑠(Ω) and with ⟨⋅,⋅⟩ the inner product in 𝐿²(Ω). Let 𝔗={𝐾} be a shape-regular quasiuniform partition of Ω of mesh size ℎ into closed elements 𝐾 [8–10]. The vector and scalar finite element spaces are 𝕎_ℎ∶=𝐰_ℎ∈𝐋²(Ω)∶𝐰_ℎ𝐾(𝐾),𝕍∈𝒫∀𝐾∈𝔗_ℎ∶=𝕎_ℎ∩𝐇¹₀(Ω),𝕋_ℎ∶=𝜃_ℎ∈𝐻¹₀(Ω)∩𝐶⁰(Ω)∶𝜃_ℎ𝐾(𝐾),ℙ∈𝒬∀𝐾∈𝔗_ℎ∶=𝑞_ℎ∈𝐿²₀(Ω)∩𝐶⁰(Ω)∶𝑞_ℎ𝐾∈ℛ(𝐾)∀𝐾∈𝔗, (3) where 𝒫(𝐾), 𝒬(𝐾), and ℛ(𝐾) are spaces of polynomials with degree bounded uniformly with respect to 𝐾∈𝔗 [9, 10]. We stress that the space ℙ_ℎ is composed of continuous functions to ensure the crucial equality ∇⋅𝐰_ℎ,𝜙_ℎ=−𝐰_ℎ,∇𝜙_ℎ,∀𝐰_ℎ∈𝕍_ℎ,∀𝜙_ℎ∈ℙ_ℎ. (4) Using the following discrete counterpart of the form 𝒩(𝐮,𝐯,𝐰)∶=⟨(𝐮⋅∇)𝐯,𝐰⟩𝒩𝐮_ℎ,𝐯_ℎ,𝐰_ℎ1∶=2𝐮_ℎ⋅∇𝐯_ℎ,𝐰_ℎ−12𝐮_ℎ⋅∇𝐰_ℎ,𝐯_ℎ, (5) we now introduce the fully discrete SGUM to solve the evolution Boussinesq equations (1).

Algorithm 1 (the fully discrete stabilized Gauge-Uzawa method). Compute 𝜃¹_ℎ, 𝐮¹_ℎ, and 𝑝¹_ℎ via any method satisfying Assumption 4 below, and set 𝜓¹_ℎ=(−2𝜏/3)𝑝¹_ℎ and 𝑞¹_ℎ=0. Repeat for 1≤𝑛≤𝑁=[𝑇/𝜏−1].

Step 1. Set 𝐮^∗_ℎ=2𝐮^𝑛_ℎ−𝐮_ℎ^𝑛−1 and 𝜃^∗_ℎ=2𝜃^𝑛_ℎ−𝜃_ℎ^𝑛−1, and then find 𝐮_ℎ^𝑛+1∈𝕍_ℎ as the solution of 13𝐮2𝜏^ℎ𝑛+1−4𝐮^𝑛_ℎ+𝐮_ℎ^𝑛−1,𝐰_ℎ+∇𝑝^𝑛_ℎ,𝐰_ℎ+𝒩𝐮^∗_ℎ,𝐮_ℎ^𝑛+1,𝐰_ℎ+𝜇∇𝐮_ℎ^𝑛+1,∇𝐰_ℎ+𝜅𝜇²𝐠𝜃^∗_ℎ,𝐰_ℎ=𝐟𝑡^𝑛+1,𝐰_ℎ,∀𝐰_ℎ∈𝕍_ℎ. (6)

Step 2. Find 𝜓_ℎ^𝑛+1∈ℙ_ℎ as the solution of ∇𝜓_ℎ^𝑛+1,∇𝜙_ℎ=∇𝜓^𝑛_ℎ,∇𝜙_ℎ+∇⋅𝐮_ℎ^𝑛+1,𝜙_ℎ,∀𝜙_ℎ∈ℙ_ℎ. (7)

Step 3. Update 𝐮_ℎ^𝑛+1 and 𝑞_ℎ^𝑛+1∈ℙ_ℎ according to 𝐮_ℎ^𝑛+1=𝐮_ℎ^𝑛+1+∇𝜓_ℎ^𝑛+1−𝜓^𝑛_ℎ𝑞_ℎ𝑛+1,,𝜙_ℎ=𝑞^𝑛_ℎ,𝜙_ℎ−∇⋅𝐮_ℎ^𝑛+1,𝜙_ℎ,∀𝜙_ℎ∈ℙ_ℎ. (8)(9)

Step 4. Update pressure 𝑝_ℎ^𝑛+1 by 𝑝_ℎ^𝑛+1=−3𝜓_ℎ^𝑛+12𝜏+𝜇𝑞_ℎ^𝑛+1. (10)

Step 5. Find 𝜃_ℎ^𝑛+1∈𝕋_ℎ as the solution of 12𝜏3𝜃_ℎ^𝑛+1−4𝜃^𝑛_ℎ+𝜃_ℎ^𝑛−1,𝜙_ℎ+𝒩𝐮_ℎ^𝑛+1,𝜃_ℎ^𝑛+1,𝜙_ℎ+𝜆𝜇∇𝜃_ℎ^𝑛+1,∇𝜙_ℎ=𝑏𝑡^𝑛+1,𝜙_ℎ,∀𝜙_ℎ∈𝕋_ℎ. (11)

Remark 2 (discontinuity of 𝐮_ℎ^𝑛+1). We note that 𝐮_ℎ^𝑛+1 is a discontinuous function across inter-element boundaries and that, in light of (7) and (8), 𝐮_ℎ^𝑛+1 is discrete divergence free in the sense that 𝐮_ℎ^𝑛+1,∇𝜙_ℎ=𝐮_ℎ^𝑛+1+𝛿∇𝜓_ℎ^𝑛+1,∇𝜙_ℎ=0,∀𝜙_ℎ∈ℙ_ℎ. (12)

We now summarize the results of this paper along with organization. We introduce appropriate Assumptions 1–5 in Section 2 and introduce well-known lemmas. In Section 3, we prove the stability result.

Theorem 3 (stability). The SGUM is unconditionally stable in the sense that, for all 𝜏>0, the following a priori bound holds: 𝐮_ℎ^𝑁+12₀+𝐮_ℎ^𝑁+12₀+2𝐮_ℎ^𝑁+1−𝐮^𝑁_ℎ²₀+2𝜃_ℎ^𝑁+1−𝜃^𝑁_ℎ²₀+𝜃_ℎ^𝑁+12₀+^𝑁_𝑛=1𝛿𝛿𝐮_ℎ^𝑛+12₀+𝛿𝛿𝜃_ℎ^𝑛+12₀+∇𝛿𝜓_ℎ^𝑛+12₀+3∇𝜓_ℎ^𝑁+12₀+2𝜇𝜏𝑞_ℎ^𝑁+12₀+𝜇𝜏^𝑁_𝑛=1∇𝐮_ℎ^𝑛+12₀+𝜆∇𝜃_ℎ^𝑛+1≤𝐮¹_ℎ²₀+2𝐮¹_ℎ−𝐮⁰_ℎ²₀+2𝜃¹_ℎ−𝜃⁰_ℎ²₀+𝜃¹_ℎ²₀+3∇𝜓¹_ℎ²₀+2𝜇𝜏𝑞¹_ℎ²₀+𝐶𝜏^𝑁_𝑛=1𝐟(𝑡^𝑛+1)²₋₁+𝑏(𝑡^𝑛+1)²₋₁. (13)

We then will carry out the following optimal error estimates through several lemmas in Section 4.

Theorem 4 (error estimates). Suppose the exact solution of (1) is smooth enough and 𝜏≤𝐶ℎ. If Assumptions 1, 3, 4, and 5 mentioned later hold, then the errors of Algorithm 1 will be bound by 𝜏^𝑁_𝑛=1𝐮(𝑡^𝑛+1)−𝐮_ℎ^𝑛+12₀+𝐮(𝑡^𝑛+1)−𝐮_ℎ^𝑛+12₀≤𝐶𝜏⁴+ℎ⁴,𝜏^𝑁_𝑛=1𝐮(𝑡^𝑛+1)−𝐮_ℎ^𝑛+12₁+ℎ²𝜃(𝑡^𝑛+1)−𝜃_ℎ^𝑛+12₁≤𝐶𝜏²+ℎ²,𝜏^𝑁_𝑛=1𝜃(𝑡^𝑛+1)−𝜃_ℎ^𝑛+12₀≤𝐶𝜏⁴+ℎ⁴. (14) Moreover, if Assumption 2 also hold, then one has 𝜏^𝑁_𝑛=1𝑝(𝑡^𝑛+1)−𝑝_ℎ^𝑛+12₀≤𝐶𝜏²+ℎ². (15)

We note that the condition 𝜏≤𝐶ℎ in Theorem 4 can be omitted for the linearized Boussinesq equations (see Remark 16). Finally, we perform numerical tests in Section 5 to check accuracy and physically relevant numerical simulations, the Bénard convection problem and the thermal driven cavity flow.

2. Preliminaries

In this section, we introduce 5 assumptions and well known lemmas to use in proof of main theorems. We resort to a duality argument for −Δ𝐯+∇𝑟=𝐳,∇⋅𝐯=0,inΩ, (16) which is the stationary Stokes system with vanishing boundary condition 𝐯=𝟎, as well as Poisson’s equation −Δ𝜔=𝜉,inΩ, (17) with boundary condition 𝜔=0.

We now state a basic assumption about Ω.

Assumption 1 (regularity of (𝐯,𝑟) and 𝜔). The unique solutions {𝐯,𝑟} of (16) and 𝜔 of (17) satisfy ‖𝐯‖₂+‖𝑟‖₁≤𝐶‖𝐳‖₀,‖𝜔‖₂≤𝐶𝜉₀. (18) We remark that the validity of Assumption 1 is known if 𝜕Ω is of class 𝐂² [11, 12], or if  𝜕Ω is a two-dimensional convex polygon [13] and is generally believed for convex polyhedral [12].

We impose the following properties for relations between the spaces 𝕍_ℎ and ℙ_ℎ.

Assumption 2 (discrete inf-sup). There exists a constant 𝛽>0 such that inf_{𝑠ℎ∈ℙℎ}sup_{𝐰ℎ∈𝕍ℎ}∇⋅𝐰_ℎ,𝑠_ℎ𝐰_ℎ1𝑠_ℎ0≥𝛽. (19)

Assumption 3 (shape regularity and quasiuniformity [8–10]). There exists a constant 𝐶>0 such that the ratio between the diameter ℎ_𝐾 of an element 𝐾∈𝔗 and the diameter of the largest ball contained in 𝐾 is bounded uniformly by 𝐶, and ℎ_𝐾 is comparable with the meshsize ℎ for all 𝐾∈𝔗.

In order to launch Algorithm 1, we need to set (𝐮¹_ℎ,𝑝¹_ℎ,𝜃¹_ℎ) via any first-order methods which hold the following conditions.

Assumption 4 (the setting of the first step values). Let (𝐮(𝑡¹),𝑝(𝑡¹),𝜃(𝑡¹)) be the exact solution of (1) at 𝑡=𝑡¹. The first step value (𝐮¹_ℎ,𝑝¹_ℎ,𝜃¹_ℎ) satisfies 𝐮𝑡¹−𝐮¹_ℎ0+𝜃𝑡¹−𝜃¹_ℎ0≤𝐶𝜏²+ℎ²,𝐮(𝑡¹)−𝐮¹_ℎ1+𝑝(𝑡¹)−𝑝¹_ℎ1+𝜃(𝑡¹)−𝜃¹_ℎ1≤𝐶(𝜏+ℎ). (20)

Assumption 5 (approximability [8–10]). For each (𝐰,𝜂,𝑠)∈𝐇²(Ω)×𝐻²(Ω)×𝐻¹(Ω), there exist approximations (𝐰_ℎ,𝜂_ℎ,𝑠_ℎ)∈𝕍_ℎ×𝕋_ℎ×ℙ_ℎ such that 𝐰−𝐰_ℎ0+ℎ𝐰−𝐰_ℎ1≤𝐶ℎ²‖𝐰‖₂,𝜂−𝜂_ℎ0+ℎ𝜂−𝜂_ℎ1≤𝐶ℎ²𝜂₂,𝑠−𝑠_ℎ0≤𝐶ℎ‖𝑠‖₁. (21)

The following elementary but crucial relations are derived in [14].

Lemma 5 (inverse inequality). If 𝐼_ℎ denotes the Clément interpolant, then 𝐼_ℎ𝐰_𝐋3(Ω)≤𝐶ℎ^−𝑑/6𝐼_ℎ𝐰₀,𝐰−𝐼_ℎ𝐰_𝐋3(Ω)≤𝐶ℎ^2−𝑑/6‖𝐰‖₂. (22)

Lemma 6 (div-grad relation). If  𝐰∈𝐇¹₀(Ω), then ‖∇⋅𝐰‖₀≤‖∇𝐰‖₀. (23)

Let now (𝐯_ℎ,𝑟_ℎ)∈𝕍_ℎ×ℙ_ℎ indicate the finite element solution of (16); namely, ∇𝐯_ℎ,∇𝐰_ℎ+∇𝑟_ℎ,𝐰_ℎ=𝐳,𝐰_ℎ,∀𝐰_ℎ∈𝕍_ℎ,∇⋅𝐯_ℎ,𝑠_ℎ=0,∀𝑠_ℎ∈ℙ_ℎ. (24)

Then we can find the well-known lemmas in [8–10].

Lemma 7 (error estimates for mixed FEM). Let (𝐯,𝑟)∈𝐇¹₀(Ω)×𝐿²₀(Ω) be the solutions of (16), and, (𝐯_ℎ,𝑟_ℎ)=𝔖_ℎ(𝐯,𝑟)∈𝕍_ℎ×ℙ_ℎ be the Stokes projection defined by (24), respectively. If Assumptions 1, 2, 3, and 5 hold, then 𝐯−𝐯_ℎ0+ℎ𝐯−𝐯_ℎ1+ℎ𝑟−𝑟_ℎ0≤𝐶ℎ²‖𝐯‖₂+‖𝑟‖₁,𝐯−𝐯_ℎ∶=𝐯−𝐯_{ℎ𝐋∞(Ω)}+∇(𝐯−𝐯_ℎ)_𝐋3(Ω)≤𝐶‖𝐳‖₀. (25)(26)

Proof. Inequality (25) is standard [8–10]. To establish (26), we just deal with the 𝐋^∞-norm, since the other can be treated similarly. If 𝐼_ℎ denotes the Clément interpolant, then ‖𝐯−𝐼_ℎ𝐯‖_𝐋^∞_(Ω)≤𝐶‖𝐯‖₂ and 𝐼_ℎ𝐯−𝐯_{ℎ𝐋∞(Ω)}≤𝐶ℎ^−𝑑/2𝐼_ℎ𝐯−𝐯_ℎ𝐋2(Ω)≤𝐶‖𝐯‖₂ (27) as a consequence of an inverse estimate. This completes the proof.

Remark 8 (𝐻¹ stability of 𝑟_ℎ). The bound ‖∇𝑟_ℎ‖₀≤𝐶(‖𝐯‖₂+‖𝑟‖₁) is a simple by-product of (16). To see this, we add and subtract 𝐼_ℎ𝑞, use the stability of 𝐼_ℎ in 𝐻¹, and observe that (25) implies ‖∇(𝑟_ℎ−𝐼_ℎ𝑟)‖₀≤𝐶ℎ⁻¹‖𝑟_ℎ−𝐼_ℎ𝑟‖₀≤𝐶.

Lemma 9 (error estimates for Poisson's equation [8, 10]). Let 𝜔∈𝐻¹₀(Ω) be the solution of Poisson's equation (17), and, 𝜔_ℎ∈𝕋_ℎ be its finite element approximation ∇𝜔_ℎ,∇𝜙_ℎ=𝜉,𝜙_ℎ,∀𝜙_ℎ∈𝕋_ℎ. (28) If Assumptions 1, 3, and 5 hold, then there exists a positive constant 𝐶 satisfying 𝜔−𝜔_ℎ0+ℎ𝜔−𝜔_ℎ1≤𝐶ℎ²‖𝜔‖₂≤𝐶ℎ²𝜉₀,𝜔−𝜔_ℎ≤𝐶𝜉₀. (29)

We finally state without proof several properties of the nonlinear form 𝒩. In view of (5), we have the following properties of 𝒩 for all 𝐮_ℎ,𝐯_ℎ,𝐰_ℎ∈𝕍_ℎ: 𝒩𝐮_ℎ,𝐯_ℎ,𝐰_ℎ=−𝒩𝐮_ℎ,𝐰_ℎ,𝐯_ℎ,𝒩𝐮_ℎ,𝐯_ℎ,𝐯_ℎ=0,∇⋅𝐮=0⇒𝒩𝐮,𝐯_ℎ,𝐰_ℎ=(𝐮⋅∇)𝐯_ℎ,𝐰_ℎ=−(𝐮⋅∇)𝐰_ℎ,𝐯_ℎ. (30)(31) Applying Sobolev imbedding lemma yields the following useful results.

Lemma 10 (bounds on nonlinear convection [1, 12]). Let 𝐮,𝐯∈𝐇²(Ω) with ∇⋅𝐮=0, and let 𝐮_h,𝐯_h,𝐰_h∈𝕍_ℎ. Then 𝒩𝐮,𝐯_ℎ,𝐰_ℎ‖𝐮‖≤𝐶¹𝐯_ℎ1𝐰_ℎ1‖𝐮‖₂𝐯_ℎ1𝐰_ℎ0‖𝐮‖₂𝐯_ℎ0𝐰_ℎ1,𝒩𝐮_ℎ,𝐯,𝐰_ℎ≤𝐶𝐮_ℎ0‖𝐯‖₂𝐰_ℎ1. (32)(33) In addition, 𝒩𝐮_ℎ,𝐯_ℎ,𝐰_ℎ𝐮≤𝐶^ℎ0𝐯_ℎ𝐰_ℎ1𝐮_ℎ𝐋3(Ω)𝐯_ℎ1𝐰_ℎ1. (34)

Remark 11 (Treatment of convection term). As we mention in Remark 2, 𝐮_ℎ^𝑛+1 is discontinuous across inter-element boundaries and so 𝐮_ℎ^𝑛+1∉𝐇¹(Ω). Thus, we can't directly apply (32) anymore to treat the convection term. To solve this difficulty, we apply (34) together with (26) and inverse inequality Lemma 5.

We will use the following algebraic identities frequently to treat time derivative terms.

Lemma 12 (inner product of time derivative terms). For any sequence {𝑧^𝑛}^𝑁_𝑛=0, one has 23𝑧^𝑛+1−4𝑧^𝑛+𝑧^𝑛−1,𝑧^𝑛+1=𝛿𝑧^𝑛+12₀+𝛿2𝑧^𝑛+1−𝑧^𝑛2₀+𝛿𝛿𝑧^𝑛+12₀,2𝑧^𝑛+1−𝑧^𝑛,𝑧ⁿ⁺¹=𝑧^𝑛+12₀−𝑧^𝑛2₀+𝑧^𝑛+1−𝑧^𝑛2₀,2𝑧^𝑛+1−𝑧^𝑛,𝑧^𝑛=𝑧^𝑛+12₀−𝑧^𝑛2₀−𝑧^𝑛+1−𝑧^𝑛2₀. (35)(36)(37)

3. Proof of Stability

We show that the SGUM is unconditionally stable via a standard energy method in this section. We start to prove stability with rewriting the momentum equation (6) by using (8) and (10) as follows: 12𝜏3𝐮_ℎ^𝑛+1−4𝐮^𝑛_ℎ+𝐮_ℎ^𝑛−1,𝐰_ℎ+𝒩𝐮^∗_ℎ,𝐮_ℎ^𝑛+1,𝐰_ℎ+𝜅𝜇²𝐠𝜃^∗_ℎ,𝐰_ℎ−∇3𝜓2𝜏^ℎ𝑛+1−𝜇𝑞^𝑛_ℎ,𝐰_ℎ+𝜇∇𝐮_ℎ^𝑛+1,∇𝐰_ℎ=𝐟𝑡^𝑛+1,𝐰_ℎ. (38) We now choose 𝐰_ℎ=4𝜏𝐮_ℎ^𝑛+1 and apply (35) to obtain 𝛿𝐮_ℎ^𝑛+12₀+𝛿2𝐮_ℎ^𝑛+1−𝐮^𝑛_ℎ²₀+𝛿𝛿𝐮_ℎ^𝑛+12₀+4𝜏𝜇∇𝐮_ℎ^𝑛+12₀=⁴_𝑖=1𝐴_𝑖, (39) where 𝐴₁∶=6∇𝜓_ℎ^𝑛+1,𝐮_ℎ^𝑛+1,𝐴₂∶=4𝜏𝜇𝑞^𝑛_ℎ,∇⋅𝐮_ℎ^𝑛+1,𝐴₃∶=4𝜏𝐟𝑡^𝑛+1,𝐮_ℎ^𝑛+1,𝐴₄∶=−4𝜏𝜅𝜇²𝐠2𝜃^𝑛_ℎ−𝜃_ℎ^𝑛−1,𝐮_ℎ^𝑛+1. (40) We give thanks to (30) for eliminating the convection term. In light of 𝐮_ℎ^𝑛+1=𝐮_ℎ^𝑛+1−∇𝛿𝜓_ℎ^𝑛+1, (12) and (36) yield 𝐴₁=−6∇𝜓_ℎ^𝑛+1,∇𝛿𝜓_ℎ^𝑛+1=−3∇𝜓_ℎ^𝑛+12₀−∇𝜓^𝑛_ℎ²₀+∇𝛿𝜓_ℎ^𝑛+12₀. (41) Before we estimate 𝐴₂, we evaluate an inequality via choosing 𝜙_ℎ=𝛿𝑞_ℎ^𝑛+1 in (9) to get 𝛿𝑞_ℎ^𝑛+12₀=−∇⋅𝐮_ℎ^𝑛+1,𝛿𝑞_ℎ^𝑛+1≤∇⋅𝐮_ℎ^𝑛+1₀𝛿𝑞_ℎ^𝑛+1₀. (42) Lemma 6 derives ‖𝛿𝑞_ℎ^𝑛+1‖²₀≤‖∇⋅𝐮_ℎ^𝑛+1‖²₀≤‖∇𝐮_ℎ^𝑛+1‖²₀, and so (9) and (37) lead us to 𝐴₂=−4𝜇𝜏𝑞^𝑛_ℎ,𝛿𝑞_ℎ^𝑛+1=−2𝜇𝜏𝑞_ℎ^𝑛+12₀−𝑞^𝑛_ℎ²₀−𝛿𝑞_ℎ^𝑛+12₀≤−2𝜇𝜏𝑞_ℎ^𝑛+12₀−𝑞^𝑛_ℎ²₀+2𝜇𝜏∇𝐮_ℎ^𝑛+12₀. (43) Clearly, we have 𝐴₃𝜏≤𝐶𝜇𝐟(𝑡^𝑛+1)²₋₁+𝜏𝜇∇𝐮_ℎ^𝑛+12₀. (44) We now attack 𝐴₄ with ‖𝐮_ℎ^𝑛+1‖²₀=‖𝐮_ℎ^𝑛+1‖²₀+‖∇𝛿𝜓_ℎ^𝑛+1‖²₀ which comes from (8) and (12). Then we arrive at 𝐴₄≤𝐶𝜅²𝜇⁴𝜏2𝜃^𝑛_ℎ−𝜃_ℎ^𝑛−12₀+𝐶𝜏𝐮_ℎ^𝑛+12₀+∇𝛿𝜓_ℎ^𝑛+12₀. (45) In conjunction with 𝒩(𝐮^∗_ℎ,𝜃_ℎ^𝑛+1,𝜃_ℎ^𝑛+1)=0 which comes from (30), if we choose 𝜙_ℎ=4𝜏𝜃_ℎ^𝑛+1 in (11) and use (35), then we obtain 𝛿𝜃_ℎ^𝑛+12₀+𝛿2𝜃_ℎ^𝑛+1−𝜃^𝑛_ℎ²₀+𝛿𝛿𝜃_ℎ^𝑛+12₀+4𝜏𝜆𝜇∇𝜃_ℎ^𝑛+1=4𝜏𝑏𝑡^𝑛+1,𝜃^𝑛_ℎ𝜏≤𝐶𝜆𝜇𝑏(𝑡^𝑛+1)²₋₁+𝜆𝜇𝜏∇𝜃^𝑛_ℎ²₀. (46) Inserting 𝐴₁−𝐴₄ back into (39), adding with (46), and then summing over 𝑛 from 1 to 𝑁 lead to (13) by help of discrete Gronwall inequality and the equality ‖𝐮_ℎ^𝑛+1‖²₀=‖𝐮_ℎ^𝑛+1‖²₀+‖∇𝛿𝜓_ℎ^𝑛+1‖²₀. So we finish the proof of Theorem 3.

4. Error Estimates

We prove Theorem 4 which is error estimates for SGUM of Algorithm 1. This proof is carried out through several lemmas. We start to prove this theorem with defining (𝐔_ℎ^𝑛+1,𝑃_ℎ^𝑛+1)∶=𝔖_ℎ(𝐮(𝑡^𝑛+1),𝑝(𝑡^𝑛+1))∈𝕍_ℎ×ℙ_ℎ to be the Stokes projection of the true solution at time 𝑡^𝑛+1. It means that (𝐔_ℎ^𝑛+1,𝑃_ℎ^𝑛+1)∈𝕍_ℎ×ℙ_ℎ is the solution of ∇𝐮𝑡^𝑛+1−𝐔_ℎ^𝑛+1,∇𝐰_ℎ+∇𝑝𝑡^𝑛+1−𝑃_ℎ^𝑛+1,𝐰_ℎ=0,∀𝐰_ℎ∈𝕍_ℎ,∇⋅𝐔_ℎ^𝑛+1,𝜙_ℎ=0,𝜙_ℎ∈ℙ_ℎ. (47) We also define Θ_ℎ^𝑛+1∈𝕋_ℎ as the solution of ∇Θ_ℎ^𝑛+1,∇𝜙_ℎ=∇𝜃𝑡^𝑛+1,∇𝜙_ℎ,∀𝜙_ℎ∈𝕋_ℎ. (48) And we denote notations 𝐆^𝑛+1∶=𝐮(𝑡^𝑛+1)−𝐔_ℎ^𝑛+1, 𝑔^𝑛+1∶=𝑝𝑡^𝑛+1−𝑃_ℎ^𝑛+1,𝜂^𝑛+1∶=𝜃𝑡^𝑛+1−Θ_ℎ^𝑛+1. (49) From Lemma 7, we can deduce 𝐆^𝑛+12₀+ℎ²𝐆^𝑛+12₁+ℎ²𝑔^𝑛+12₀≤𝐶ℎ⁴𝐮(𝑡^𝑛+1)²₂+𝑝(𝑡^𝑛+1)²₁,𝛿𝐆^𝑛+12₀+ℎ²𝛿𝐆^𝑛+12₁+ℎ²𝛿𝑔^𝑛+12₀≤𝐶𝜏ℎ^4𝑡𝑛+1_𝑡𝑛𝐮_𝑡(𝑡)²₂+𝑝_𝑡(𝑡)²₁𝑑𝑡. (50)(51) In conjunction with the definition of |‖⋅‖| in (26), we can derive 𝐆^𝑛+1≤𝐶. (52) We now carry out error evaluation by comparing (65) with (6)–(10) and then by comparing (82) and (11). We derive strong estimates of order 1, and this result is instrumental in proving weak estimates of order 2 for the errors 𝐄_ℎ^𝑛+1∶=𝐔_ℎ^𝑛+1−𝐮_ℎ^𝑛+1,𝐄_ℎ^𝑛+1∶=𝐔_ℎ^𝑛+1−𝐮_ℎ^𝑛+1,𝑒_ℎ^𝑛+1∶=𝑃_ℎ^𝑛+1−𝑝_ℎ^𝑛+1,𝜗_ℎ^𝑛+1∶=Θ_ℎ^𝑛+1−𝜃_ℎ^𝑛+1. (53) Then, in conjunction with (50), we can readily get the same accuracy for the errors 𝐄^𝑛+1∶=𝐮𝑡^𝑛+1−𝐮_ℎ^𝑛+1,𝐄^𝑛+1∶=𝐮𝑡^𝑛+1−𝐮_ℎ^𝑛+1,𝑒^𝑛+1∶=𝑝𝑡^𝑛+1−𝑝_ℎ^𝑛+1,𝜗^𝑛+1∶=𝜃𝑡^𝑛+1−𝜃_ℎ^𝑛+1. (54) Additionally, we denote 𝜀_ℎ^𝑛+1∶=𝑃_ℎ^𝑛+1+3𝜓_ℎ^𝑛+1.2𝜏 (55) We readily obtain the following properties: 𝐆^𝑛+1=𝐄^𝑛+1−𝐄_ℎ^𝑛+1=𝐄^𝑛+1−𝐄_ℎ^𝑛+1,𝐄_ℎ^𝑛+1=𝐄_ℎ^𝑛+1+∇𝛿𝜓_ℎ^𝑛+1,𝐄^𝑛+1=𝐄_ℎ^𝑛+1=𝐆^𝑛+1=𝟎,on𝜕Ω. (56) Moreover, from (12), 𝐄_ℎ^𝑛+1,∇𝜙_ℎ=𝐄^𝑛+1,∇𝜙_ℎ=𝐆^𝑛+1,∇𝜙_ℎ=0,∀𝜙_ℎ∈ℙ_ℎ. (57) Whence we deduce crucial orthogonality properties: 𝐄_ℎ^𝑛+12₀=𝐄_ℎ^𝑛+12₀+∇𝛿𝜓_ℎ^𝑛+12₀. (58) We also point out that, owing to Lemma 6, 𝑞_ℎ^𝑛+1∈ℙ_ℎ defined in (9) satisfies 𝑞_ℎ^𝑛+1−𝑞^𝑛_ℎ0≤∇𝐄^𝑛+1₀. (59) In conjunction with 𝜂^𝑛+1=𝜃(𝑡^𝑛+1)−Θ_ℎ^𝑛+1=𝜗^𝑛+1−𝜗_ℎ^𝑛+1, Assumption 5 leads to 𝜂^𝑛+12₀+ℎ²𝜂^𝑛+1₁≤𝐶ℎ⁴𝜃(𝑡^𝑛+1)²₂,𝛿𝜂^𝑛+12₀+ℎ²𝛿𝜂^𝑛+12₁≤𝐶𝜏ℎ^4𝑡𝑛+1_𝑡𝑛𝜃_𝑡2𝑑𝑡. (60)

We now estimate the first-order accuracy for velocity and temperature in Lemma 13 and then the 2nd-order accuracy for time derivative of velocity and temperature in Lemma 15. The result of Lemma 13 is instrumental to treat the convection term in proof of Lemma 15. We will use Lemmas 13 and 15 to prove optimal error accuracy in Lemma 17. Finally, we will prove pressure error estimate in Lemma 18.

Lemma 13 (reduced rate of convergence for velocity and temperature). Suppose the exact solution of (1) is smooth enough. If Assumptions 3–5 hold, then the velocity and temperature error functions satisfy 𝐄_ℎ^𝑁+12₀+𝐄_ℎ^𝑁+12₀+2𝐄_ℎ^𝑁+1−𝐄^𝑁_ℎ²₀+2𝜇𝜏𝑞_ℎ^𝑁+12₀+4𝜏²3∇𝜀_ℎ^𝑁+12₀+𝜗_ℎ^𝑁+12₀+2𝜗_ℎ^𝑁+1−𝜗^𝑁_ℎ²₀+𝜇𝜏^𝑁_𝑛=1∇𝐄_ℎ^𝑛+12₀+𝜆∇𝜗_ℎ^𝑛+12₀+^𝑁_𝑛=1𝛿𝛿𝐄_ℎ^𝑛+12₀+∇𝛿𝜓_ℎ^𝑛+12₀+𝛿𝛿𝜗_ℎ^𝑛+12₀≤𝐶𝜏²+ℎ². (61)

Proof. We resort to the Taylor theorem to write (1) as follows: 3𝐮𝑡^𝑛+1−4𝐮𝑡^𝑛+𝐮𝑡^𝑛−1+𝐮𝑡2𝜏_𝑛+1⋅∇𝐮𝑡^𝑛+1+∇𝑝𝑡^𝑛+1−𝜇Δ𝐮𝑡^𝑛+1+𝜅𝜇²𝐠𝜃𝑡^𝑛+1=𝐑^𝑛+1+𝐟𝑡^𝑛+1,3𝜃𝑡^𝑛+1−4𝜃𝑡^𝑛+𝜃𝑡^𝑛−1𝑡2𝜏+𝐮^𝑛+1⋅∇𝜃𝑡^𝑛+1−𝜆𝜇Δ𝜃𝑡^𝑛+1=𝑄^𝑛+1+𝑏𝑡^𝑛+1, (62)(63) where 𝐑^𝑛+11∶=𝜏^𝑡𝑛+1_𝑡^𝑛𝐮_𝑡𝑡𝑡(𝑠)𝑠−𝑡^𝑛2−1𝑑𝑠4𝜏^𝑡𝑛+1_𝑡^𝑛−1𝐮_𝑡𝑡𝑡(𝑠)𝑡^𝑛−1−𝑠²𝑄𝑑𝑠,_𝑛+11∶=𝜏^𝑡𝑛+1_𝑡^𝑛𝜃_𝑡𝑡𝑡(𝑠)𝑠−𝑡^𝑛2−1𝑑𝑠4𝜏^𝑡𝑛+1_𝑡^𝑛−1𝜃_𝑡𝑡𝑡(𝑠)𝑡^𝑛−1−𝑠²𝑑𝑠 (64) are the truncation errors. In conjunction with the definition of the Stokes projection {𝐔_ℎ^𝑛+1,𝑃_ℎ^𝑛+1}, we readily get a weak formulation for (62), ∀𝐰_ℎ∈𝕍_ℎ: 1𝑡2𝜏3𝐮^𝑛+1−4𝐮𝑡^𝑛+𝐮𝑡^𝑛−1,𝐰_ℎ+𝒩𝐮𝑡^𝑛+1,𝐮𝑡^𝑛+1,𝐰_ℎ+∇𝑃_ℎ^𝑛+1,𝐰_ℎ+𝜇∇𝐔_ℎ^𝑛+1,∇𝐰_ℎ+𝜅𝜇²𝐠𝜃𝑡^𝑛+1,𝐰_ℎ=𝐑^𝑛+1,𝐰_ℎ+𝐟𝑡^𝑛+1,𝐰_ℎ. (65) We replace the pressure 𝑝^𝑛_ℎ term in (6) with (10) and then subtract it from (65) to obtain 13𝐄2𝜏^𝑛+1−4𝐄^𝑛+𝐄^𝑛−1,𝐰_ℎ+𝜇𝜏∇𝐄_ℎ^𝑛+1,∇𝐰_ℎ=−∇𝛿𝑃_ℎ^𝑛+1+𝜀^𝑛_ℎ,𝐰_ℎ+𝜇∇𝑞^𝑛_ℎ,𝐰_ℎ+𝒩2𝐮^𝑛_ℎ−𝐮_ℎ^𝑛−1,𝐮_ℎ^𝑛+1,𝐰_ℎ−𝒩𝐮𝑡^𝑛+1,𝐮𝑡^𝑛+1,𝐰_ℎ−𝜅𝜇²𝐠𝜃𝑡^𝑛+1−2𝜃^𝑛_ℎ+𝜃_ℎ^𝑛−1,𝐰_ℎ+𝐑^𝑛+1,𝐰_ℎ. (66) Choosing 𝐰_ℎ=4𝜏𝐄_ℎ^𝑛+1=4𝜏(𝐄^𝑛+1−𝐆^𝑛+1)=4𝜏(𝐄_ℎ^𝑛+1+∇𝛿𝜓_ℎ^𝑛+1) in (66) and using (58) and (35) yield 𝛿𝐄_ℎ^𝑛+12₀+𝛿2𝐄_ℎ^𝑛+1−𝐄^𝑛_ℎ²₀+𝛿𝛿𝐄_ℎ^𝑛+12₀+6∇𝛿𝜓_ℎ^𝑛+12₀+4𝜇𝜏∇𝐄_ℎ^𝑛+12₀=⁷_𝑛=1𝐴_𝑖, (67) where 𝐴₁∶=−4𝜏𝒩𝐮𝑡^𝑛+1,𝐮𝑡^𝑛+1,𝐄_ℎ^𝑛+1+4𝜏𝒩2𝐮^𝑛_ℎ−𝐮_ℎ^𝑛−1,𝐮_ℎ^𝑛+1,𝐄_ℎ^𝑛+1,𝐴₂∶=−23𝐆^𝑛+1−4𝐆^𝑛+𝐆^𝑛−1,𝐄_ℎ^𝑛+1,𝐴₃∶=−4𝜏∇𝛿𝑃_ℎ^𝑛+1,𝐄_ℎ^𝑛+1,𝐴₄∶=−4𝜏∇𝜀^𝑛_ℎ,𝐄_ℎ^𝑛+1,𝐴₅∶=4𝜇𝜏∇𝑞^𝑛_ℎ,𝐄_ℎ^𝑛+1,𝐴₆∶=4𝜏𝐑^𝑛+1,𝐄_ℎ^𝑛+1,𝐴₇∶=−4𝜏𝜅𝜇²𝐠𝜃𝑡^𝑛+1−2𝜃^𝑛_ℎ+𝜃_ℎ^𝑛−1,𝐄_ℎ^𝑛+1. (68) We now estimate all the terms from 𝐴₁ to 𝐴₇, respectively.  To tackle 𝐴₁, we first add and subtract 2𝐮(𝑡^𝑛)−𝐮(𝑡^𝑛−1) to obtain 𝐴₁=−4𝜏𝒩𝛿𝛿𝐮𝑡^𝑛+1,𝐮𝑡^𝑛+1,𝐄_ℎ^𝑛+1−4𝜏𝒩2𝐮^𝑛_ℎ−𝐮_ℎ^𝑛−1,𝐄^𝑛+1,𝐄_ℎ^𝑛+1−4𝜏𝒩2𝐄^𝑛−𝐄^𝑛−1,𝐮𝑡^𝑛+1,𝐄_ℎ^𝑛+1. (69) Because of 𝒩(2𝐮^𝑛_ℎ−𝐮_ℎ^𝑛−1,𝐄_ℎ^𝑛+1,𝐄_ℎ^𝑛+1)=0, which comes from (30), the second term of 𝐴₁ can be replaced by 4𝜏𝒩2𝐄^𝑛−𝐄^𝑛−1−2𝐮𝑡^𝑛+𝐮𝑡^𝑛−1,𝐆^𝑛+1,𝐄_ℎ^𝑛+1. (70) If we apply Lemma 10, then we can readily obtain 𝐴₁≤𝐶𝜏𝛿𝛿𝐮𝑡^𝑛+1₀𝐮𝑡^𝑛+1₂+2𝐄^𝑛−𝐄^𝑛−1₀𝐮𝑡^𝑛+1₂𝐄_ℎ^𝑛+1₁+𝐶𝜏2𝐄^𝑛−𝐄^𝑛−1₀𝐆^𝑛+1+2𝐮(𝑡^𝑛)−𝐮(𝑡^𝑛−1)₂𝐆^𝑛+1₀𝐄_ℎ^𝑛+1₁. (71) Since we have ‖𝐮(𝑡^𝑛+1)‖₂+‖|𝐆^𝑛+1|‖≤𝑀 according to (52), we arrive at 𝐴₁≤𝐶𝜏2𝐄^𝑛_ℎ−𝐄_ℎ^𝑛−12₀+2𝐆^𝑛−𝐆^𝑛−12₀+G^𝑛+12₀+𝜇𝜏∇𝐄_ℎ^𝑛+12₀+𝐶𝜏⁴𝜇×^𝑡𝑛+1_𝑡𝑛−1𝐮_𝑡𝑡𝑡(𝑡)²₀𝑑𝑡. (72) In light of (51) and (58), 𝐴₂ becomes 𝐴₂≤𝐶𝜏𝐄_ℎ^𝑛+12₀+14∇𝛿𝜓_ℎ^𝑛+12₀+𝐶ℎ^4𝑡𝑛+1_𝑡𝑛−1𝐮_𝑡(𝑡)²₂+𝑝_𝑡(𝑡)²₁𝑑𝑡. (73) In order to estimate 𝐴₃ and 𝐴₄, we note that the cancellation law (12) gives ⟨∇𝜀^𝑛_ℎ,𝐄_ℎ^𝑛+1⟩=0. Then 𝐄_ℎ^𝑛+1=𝐄_ℎ^𝑛+1+∇𝛿𝜓_ℎ^𝑛+1 yields 𝐴₃=−4𝜏∇𝛿𝑃_ℎ^𝑛+1,∇𝛿𝜓_ℎ^𝑛+1≤14∇𝛿𝜓_ℎ^𝑛+12₀+𝐶𝜏³×^𝑡𝑛+1_𝑡𝑛𝐮_𝑡(𝑡)²₂+𝑝_𝑡(𝑡)²₁𝑑𝑡. (74) In conjunction with the definition 𝜀_ℎ^𝑛+1=𝑃_ℎ^𝑛+1+3𝜓_ℎ^𝑛+1/2𝜏, 𝐴₄ can be evaluated by 𝐴₄=−4𝜏∇𝜀^𝑛_ℎ,∇𝛿𝜓_ℎ^𝑛+1=−8𝜏²3∇𝜀^𝑛_ℎ,∇𝛿𝜀_ℎ^𝑛+1−𝛿𝑃_ℎ^𝑛+1≤−4𝜏²3∇𝜀_ℎ^𝑛+12₀−∇𝜀^𝑛_ℎ²₀−∇𝛿𝜀_ℎ^𝑛+12₀+𝐶𝜏³∇𝜀^𝑛_ℎ²₀+𝐶𝜏∇𝛿𝑃_ℎ^𝑛+12₀. (75) If we now apply inequality (𝑎+𝑏)²≤4𝑎²+(4/3)𝑏², then we can get 4𝜏²3∇𝛿𝜀_ℎ^𝑛+12₀=4𝜏²3∇𝛿𝑃_ℎ^𝑛+1+32𝜏∇𝛿𝜓_ℎ^𝑛+12₀≤𝐶𝜏²∇𝛿𝑃_ℎ^𝑛+12₀+4∇𝛿𝜓_ℎ^𝑛+12₀. (76) So we arrive at 𝐴₄≤−4𝜏²3∇𝜀_ℎ^𝑛+12₀−∇𝜀^𝑛_ℎ²₀+𝐶𝜏³∇𝜀^𝑛_ℎ²₀+4∇𝛿𝜓_ℎ^𝑛+12₀+𝐶𝜏^2𝑡𝑛+1_𝑡𝑛𝐮_𝑡(𝑡)²₂+𝑝_𝑡(𝑡)²₁𝑑𝑡. (77) In light of ∇⋅𝐮(𝑡^𝑛+1)=0 and (37), (9) and (59) yield 𝐴₅=4𝜇𝜏𝑞^𝑛_ℎ,∇⋅𝐮_ℎ^𝑛+1=−4𝜇𝜏𝑞^𝑛_ℎ,𝛿𝑞_ℎ^𝑛+1=−2𝜇𝜏𝑞_ℎ^𝑛+12₀−𝑞^𝑛_ℎ²₀−𝛿𝑞_ℎ^𝑛+12₀≤−2𝜇𝜏𝑞_ℎ^𝑛+12₀−𝑞^𝑛_ℎ²₀+2𝜇𝜏∇𝐄_ℎ^𝑛+12₀. (78) Also we readily get 𝐴₆≤𝐶𝜏𝐄_ℎ^𝑛+12₀+14∇𝛿𝜓_ℎ^𝑛+12₀+𝐶𝜏^4𝑡𝑛+1_𝑡𝑛−1𝐮_𝑡𝑡𝑡(𝑡)²₀𝑑𝑡. (79) The Hölder inequality and (60) yield 𝐴₇=−4𝜅𝜇²𝜏𝐠𝛿𝛿𝜃𝑡^𝑛+1−2𝜗^𝑛+𝜗^𝑛−1,𝐄_ℎ^𝑛+1≤𝐶𝜅²𝜇⁴𝜏2𝜗^𝑛_ℎ−𝜗_ℎ^𝑛−12₀+2𝜂^𝑛−𝜂^𝑛−12₀+𝜏𝐄_ℎ^𝑛+12₀+14∇𝛿𝜓_ℎ^𝑛+12₀+𝐶𝜏^4𝑡𝑛_𝑡⁺¹_𝑛−1𝜃_𝑡𝑡(𝑡)²₀𝑑𝑡. (80) Inserting the above estimates into (67) gives 𝛿𝐄_ℎ^𝑛+12₀+𝛿2𝐄_ℎ^𝑛+1−𝐄^𝑛_ℎ²₀+𝛿𝛿𝐄_ℎ^𝑛+12₀+∇𝛿𝜓_ℎ^𝑛+12₀+𝜇𝜏∇𝐄_ℎ^𝑛+12₀+4𝜏²3∇𝜀_ℎ^𝑛+12₀+2𝜇𝜏𝑞_ℎ^𝑛+12₀≤4𝜏²3∇𝜀^𝑛_ℎ²₀+𝐶𝜏³∇𝜀^𝑛_ℎ²₀+2𝜇𝜏𝑞^𝑛_ℎ²₀+𝐶𝜏𝐄_ℎ^𝑛+12₀+2𝐄^𝑛_ℎ−𝐄_ℎ^𝑛−12₀+2𝐆^𝑛−𝐆^𝑛−12₀+𝐆^𝑛+12₀+𝐶𝜅²𝜇⁴𝜏2𝜂^𝑛−𝜂^𝑛−12₀+𝐶𝜏^4𝑡𝑛_𝑡⁺¹_𝑛−1𝐮_𝑡𝑡𝑡(𝑡)²₀+𝜃_𝑡𝑡(𝑡)²₀𝑑𝑡+𝐶𝜅²𝜇⁴𝜏2𝜗^𝑛_ℎ−𝜗_ℎ^𝑛−12₀+𝐶𝜏²+ℎ^4𝑡𝑛+1_𝑡𝑛−1𝐮_𝑡(𝑡)²₂+𝑝_𝑡(𝑡)²₁𝑑𝑡. (81) On the other hand, the definition (48) of Θ_ℎ^𝑛+1∈𝕋_ℎ leads to a weak formula of (63) as, for all 𝜙_ℎ∈𝕋_ℎ, 1𝑡2𝜏3𝜃^𝑛+1−4𝜃𝑡^𝑛+𝜃𝑡^𝑛−1,𝜙_ℎ+𝒩𝐮𝑡^𝑛+1,𝜃𝑡^𝑛+1,𝜙_ℎ+𝜆𝜇∇Θ_ℎ^𝑛+1,∇𝜙_ℎ=𝑄^𝑛+1+𝑏𝑡^𝑛+1,𝜙_ℎ. (82) We now subtract (11) from (82) to derive 12𝜏3𝜗^𝑛+1−4𝜗^𝑛+𝜗^𝑛−1,𝜙_ℎ+𝒩𝐮𝑡^𝑛+1,𝜃𝑡^𝑛+1,𝜙_ℎ−𝒩𝐮_ℎ^𝑛+1,𝜃_ℎ^𝑛+1,𝜙_ℎ+𝜆𝜇∇𝜗_ℎ^𝑛+1,∇𝜙_ℎ=𝑄^𝑛+1,𝜙_ℎ. (83) Choosing 𝜙_ℎ=4𝜏𝜗_ℎ^𝑛+1=4𝜏(𝜗^𝑛+1−𝜂^𝑛+1) yields 𝛿𝜗_ℎ^𝑛+12₀+𝛿2𝜗_ℎ^𝑛+1−𝜗^𝑛_ℎ²₀+𝛿𝛿𝜗_ℎ^𝑛+12₀+4𝜆𝜇𝜏∇𝜗_ℎ^𝑛+12₀=𝐴₈+𝐴₉, (84) where 𝐴₈∶=−4𝜏𝒩𝐮𝑡^𝑛+1,𝜃𝑡^𝑛+1,𝜗_ℎ^𝑛+1+4𝜏𝒩𝐮_ℎ^𝑛+1,𝜃_ℎ^𝑛+1,𝜗_ℎ^𝑛+1,𝐴₉∶=−23𝜂^𝑛+1−4𝜂^𝑛+𝜂^𝑛−1,𝜗_ℎ^𝑛+1+4𝜏𝑄^𝑛+1,𝜗_ℎ^𝑛+1. (85) To estimate 𝐴₈, we note 𝒩(𝐄^𝑛+1−𝐮(𝑡^𝑛+1),𝜗_ℎ^𝑛+1,𝜗_ℎ^𝑛+1)=0 which comes from (30), then ‖𝐄^𝑛+1‖_𝐋³_(Ω)≤𝐶ℎ^−1/2‖𝐄^𝑛+1‖₀ and ‖𝜂^𝑛+1‖₁≤𝐶ℎ‖𝜃(𝑡^𝑛+1)‖₂ yield 𝐴₈=−4𝜏𝒩𝐄^𝑛+1,𝜃𝑡^𝑛+1,𝜗_ℎ^𝑛+1+4𝜏𝒩𝐄^𝑛+1−𝐮𝑡^𝑛+1,𝜂^𝑛+1,𝜗_ℎ^𝑛+1≤𝐶𝜏𝐄^𝑛+1₀𝜃(𝑡^𝑛+1)₂𝜗_ℎ^𝑛+1₁+𝐶𝜏𝐄^𝑛+1_𝐋3(Ω)𝜂^𝑛+1₁+𝐮(𝑡^𝑛+1)₂𝜂^𝑛+1₀𝜗_ℎ^𝑛+1₁≤𝐶𝜏𝐄𝜆𝜇^𝑛+120+𝜂^𝑛+12₀+𝜆𝜇𝜏∇𝜗_ℎ^𝑛+12₀. (86) We use ‖𝛿𝜂^𝑛+1‖²₀≤𝐶𝜏ℎ⁴∫^𝑡𝑛+1_𝑡𝑛‖𝜃_𝑡(𝑡)‖²₂𝑑𝑡 which is (60) to attack 𝐴₉: 𝐴₉≤𝐶3𝜂^𝑛+1−4𝜂^𝑛+𝜂^𝑛−1₀𝜗_ℎ^𝑛+1₀+𝐶𝜏𝑄^𝑛+1₀𝜗_ℎ^𝑛+1₀≤𝐶𝜏𝜗_ℎ^𝑛+12₀+𝐶ℎ^4𝑡𝑛+1_𝑡𝑛−1𝜃_𝑡(𝑡)²₂𝑑𝑡+𝐶𝜏^4𝑡𝑛+1_𝑡𝑛−1𝜃_𝑡𝑡𝑡(𝑡)²₀𝑑𝑡. (87) Inserting the above estimates into (84) gives 𝛿𝜗_ℎ^𝑛+12₀+𝛿2𝜗_ℎ^𝑛+1−𝜗^𝑛_ℎ²₀+𝛿𝛿𝜗_ℎ^𝑛+12₀+3𝜆𝜇𝜏∇𝜗_ℎ^𝑛+12₀≤𝐶𝜏𝜗_ℎ^𝑛+12₀+𝐶𝜏𝐄𝜆𝜇^𝑛+120+𝜂^𝑛+12₀+𝐶ℎ⁴+𝜏^4𝑡𝑛+1_𝑡𝑛−1𝜃_𝑡𝑡𝑡(𝑡)²₀+𝜃_𝑡(𝑡)²₂𝑑𝑡. (88) Adding 2 inequalities (81) and (88) and summing over 𝑛 from 1 to 𝑁 imply 𝐄_ℎ^𝑁+12₀+2𝐄_ℎ^𝑁+1−𝐄^𝑁_ℎ²₀+𝜗_ℎ^𝑁+12₀+2𝜗_ℎ^𝑁+1−𝜗^𝑁_ℎ²₀+2𝜇𝜏𝑞_ℎ^𝑛+12₀+4𝜏²3∇𝜀_ℎ^𝑛+12₀+^𝑁_𝑛=1𝛿𝛿𝐄_ℎ^𝑁+12₀+𝛿𝛿𝜗_ℎ^𝑛+12₀+∇𝛿𝜓_ℎ^𝑛+12₀+𝜇𝜏^𝑁_𝑛=1∇𝐄_ℎ^𝑛+12₀+3𝜆∇𝜗_ℎ^𝑛+12₀≤𝐄¹_ℎ²₀+2𝐄¹_ℎ−𝐄⁰_ℎ²₀+𝜗¹_ℎ²₀+2𝜗¹_ℎ−𝜗⁰_ℎ²₀+4𝜏²3∇𝜀¹_ℎ²₀+2𝜇𝜏𝑞¹_ℎ²₀+𝜏^𝑁_𝑛=1𝐄^𝑛+12₀+𝜂^𝑛+12₀+𝐶𝜏^𝑁_𝑛=1𝜗_ℎ^𝑛+12₀+2𝜗^𝑛_ℎ−𝜗_ℎ^𝑛−12₀+2𝜂^𝑛−𝜂^𝑛−12₀+𝐶𝜏³∇𝜀^𝑛_ℎ²₀+𝐶ℎ⁴+𝜏⁴×^𝑡𝑛+1_𝑡𝑛−1𝐮_𝑡𝑡𝑡(𝑡)²₀+𝜃_𝑡𝑡𝑡(𝑡)²₀+𝜃_𝑡(𝑡)²₂𝜏𝑑𝑡+𝐶²+ℎ^4𝑡𝑛+1_𝑡𝑛−1𝐮_𝑡(𝑡)²₂+𝑝_𝑡(𝑡)²₁𝑑𝑡. (89) Because of 𝑞¹=0, we have 𝜀¹_ℎ=𝑃¹_ℎ+(3/2𝜏)𝜓¹_ℎ=𝑒¹_ℎ, and thus, Assumption 4 yields ‖∇𝜀¹_ℎ‖²₀≤𝐶 and the first 4 terms in the right-hand side can be bound by Assumption 4 and properties 𝐄⁰=𝟎 and 𝑞¹=0 which are directly deduced from the conditions in Algorithm 1. Moreover, the remaining terms can be treated by the discrete Gronwall lemma. Finally, in conjunction with (58), we conclude the desired result and complete this proof.

Remark 14 (optimal estimation). In order to get optimal accuracy, we must get rid of the terms of 𝐴₃ and 𝐴₄ by applying duality argument in Lemma 17. To do this, we first evaluate the errors for time derivative of velocity and temperature in Lemma 15. Thus, we need to evaluate optimal initial errors for the case 𝑛=1, and so we have to recompute again (77). We start to rewrite 𝐴₄ as 𝐴₄=−8𝜏²3∇𝜀¹_ℎ,∇𝛿𝜀²_ℎ−𝛿𝑃²_ℎ≤−4𝜏²3∇𝜀²_ℎ²₀−∇𝜀¹_ℎ²₀−∇𝛿𝜀²_ℎ²₀+𝐶𝜏²∇𝜀¹_ℎ²₀+𝐶𝜏²∇𝛿𝑃_ℎ^𝑛+12₀≤−4𝜏²3∇𝜀²_ℎ²₀+𝐶𝜏²∇𝜀¹_ℎ²₀+4∇𝛿𝜓²_ℎ²₀+𝐶𝜏³×^𝑡2_𝑡0𝑝_𝑡𝑡(𝑡)²₁+𝐮_𝑡𝑡(𝑡)²₂𝑑𝑡. (90) In light of Assumption 4, we arrive at 𝐄²_ℎ²₀+𝐄²_ℎ²₀+2𝐄²_ℎ−𝐄¹_ℎ²₀+12𝛿𝛿𝐄²_ℎ²₀+∇𝛿𝜓²_ℎ²₀+2𝜗²_ℎ−𝜗²_ℎ²₀+𝜗²_ℎ²₀+𝜇𝜏∇𝐄²_ℎ²₀+2𝜇𝜏𝑞²_ℎ²₀+4𝜏²3∇𝜀²_ℎ²₀+𝜇𝜆𝜏∇𝜗²_ℎ²₀≤𝐶𝜏⁴+ℎ⁴. (91)

We now start to estimate errors for time derivative of velocity.

Lemma 15 (error estimate for time derivative of velocity and temperature). Suppose the exact solution of (1) is smooth enough and 𝜏≤𝐶ℎ. If Assumptions 3–5 hold, then the time derivative velocity and temperature error functions satisfy 𝛿𝐄_ℎ^𝑁+12₀+𝛿𝐄_ℎ^𝑁+12₀+2𝛿𝐄_ℎ^𝑁+1−𝛿𝐄^𝑁_ℎ²₀+2𝜇𝜏𝛿𝑞_ℎ^𝑁+12₀+𝛿𝜗_ℎ^𝑁+12₀+2𝛿𝜗_ℎ^𝑁+1−𝛿𝜗^𝑁_ℎ²₀+^𝑁_𝑛=2𝛿𝛿𝛿𝐄_ℎ^𝑛+12₀+∇𝛿𝛿𝜓_ℎ^𝑛+12₀+𝛿𝛿𝛿𝜗_ℎ^𝑛+12₀+4𝜏²3∇𝛿𝜀_ℎ^𝑁+12₀+𝜇𝜏^𝑁_𝑛=2∇𝛿𝐄_ℎ^𝑛+12₀+𝜆∇𝛿𝜗_ℎ^𝑛+12₀≤𝐶𝜏²+ℎ²². (92)

Remark 16 (the condition 𝜏≤𝐶ℎ). The assumption 𝜏≤𝐶ℎ requires to control convection terms which are used at only (100) and (117), so we can omit this condition for the linearized Boussinesq equations. However, this condition cannot be removed for nonlinear equation case, because (99) must be bounded by ℎ.

Proof. We subtract two consecutive formulas of (66) and impose 𝐰_ℎ=4𝜏𝛿𝐄_ℎ^𝑛+1 to obtain 𝛿𝐄_ℎ^𝑛+12₀+2𝛿𝐄_ℎ^𝑛+1−𝛿𝐄^𝑛_ℎ²₀+𝛿𝛿𝛿𝐄_ℎ^𝑛+12₀+6∇𝛿𝛿𝜓_ℎ^𝑛+12₀−𝛿𝐄^𝑛_ℎ²₀−2𝛿𝐄^𝑛_ℎ−𝛿𝐄_ℎ^𝑛−12₀+4𝜇𝜏∇𝛿𝐄_ℎ^𝑛+12₀=⁷_𝑛=1𝐴_𝑖, (93) where 𝐴₁∶=4𝜏𝒩𝐮𝑡^𝑛,𝐮𝑡^𝑛,𝛿𝐄_ℎ^𝑛+1−4𝜏𝒩2𝐮_ℎ^𝑛−1−𝐮_ℎ^𝑛−2,𝐮^𝑛_ℎ,𝛿𝐄_ℎ^𝑛+1−4𝜏𝒩𝐮𝑡^𝑛+1,𝐮𝑡^𝑛+1,𝛿𝐄_ℎ^𝑛+1+4𝜏𝒩2𝐮^𝑛_ℎ−𝐮_ℎ^𝑛−1,𝐮_ℎ^𝑛+1,𝛿𝐄_ℎ^𝑛+1,𝐴₂∶=−23𝛿𝐆^𝑛+1−4𝛿𝐆^𝑛+𝛿𝐆^𝑛−1,𝛿𝐄_ℎ^𝑛+1,𝐴₃∶=−4𝜏∇𝛿𝛿𝑃_ℎ^𝑛+1,𝛿𝐄_ℎ^𝑛+1,𝐴₄∶=−4𝜏∇𝛿𝜀^𝑛_ℎ,𝛿𝐄_ℎ^𝑛+1,𝐴₅∶=4𝜇𝜏∇𝛿𝑞^𝑛_ℎ,𝛿𝐄_ℎ^𝑛+1,𝐴₆∶=4𝜏𝛿𝐑^𝑛+1,𝛿𝐄_ℎ^𝑛+1,𝐴₇∶=4𝜅𝜇²𝜏×𝐠𝛿𝜃𝑡^𝑛+1−2𝜃^𝑛+𝜃^𝑛−1,𝛿𝐄_ℎ^𝑛+1. (94) We now estimate each term from 𝐴₁ to 𝐴₇  separately. The convection term 𝐴₁ can be rewritten as follows: 𝐴₁=4𝜏𝒩𝛿𝛿𝐮𝑡^𝑛,𝐮𝑡^𝑛,𝛿𝐄^𝑛+1−4𝜏𝒩𝛿𝛿𝐮𝑡^𝑛+1,𝐮𝑡^𝑛+1,𝛿𝐄^𝑛+1+4𝜏𝒩2𝐄^𝑛−1−𝐄^𝑛−2,𝐮𝑡^𝑛,𝛿𝐄^𝑛+1−4𝜏𝒩2𝐄^𝑛−𝐄^𝑛−1,𝐮𝑡^𝑛+1,𝛿𝐄^𝑛+1+4𝜏𝒩2𝐮_ℎ^𝑛−1−𝐮_ℎ^𝑛−2,𝐄^𝑛,𝛿𝐄^𝑛+1−4𝜏𝒩2𝐮^𝑛_ℎ−𝐮_ℎ^𝑛−1,𝐄^𝑛+1,𝛿𝐄^𝑛+1. (95) In estimating convection terms, we will use Lemma 10 frequently without notice. We recall ‖𝐮(𝑡)‖₂≤𝐶 to obtain 𝐴_1,1+𝐴_1,2≤𝐶𝜏𝛿𝛿𝐮(𝑡^𝑛+1)₀𝐮(𝑡^𝑛+1)₂+𝛿𝛿𝐮(𝑡^𝑛)₀𝐮(𝑡^𝑛)₂𝛿𝐄_ℎ^𝑛+1₁≤𝜇𝜏6∇𝛿𝐄_ℎ^𝑛+12₀+𝐶𝜏⁴𝜇^𝑡𝑛+1_𝑡^𝑛−2𝐮_𝑡𝑡(𝑡)²₀𝑑𝑡. (96) The result in Lemma 13, ‖2𝐄^𝑛−𝐄^𝑛−1‖₀≤𝐶(𝜏+ℎ), is essential to treat the next 2 convection terms. Invoking (33), we have 𝐴_1,3+𝐴_1,4≤𝐶𝜏2𝐄^𝑛−𝐄^𝑛−1₀×𝛿𝐮(𝑡^𝑛+1)₂𝛿𝐄_ℎ^𝑛+1₁+𝐶𝜏2𝛿𝐄^𝑛−𝛿𝐄^𝑛−1₀×𝐮(𝑡^𝑛)₂𝛿𝐄_ℎ^𝑛+1₁≤𝜇𝜏6∇𝛿𝐄_ℎ^𝑛+12₀+𝐶𝜏𝜇2𝛿𝐄^𝑛−𝛿𝐄^𝑛−12₀+𝐶𝜏²𝜏²+ℎ²𝜇×^𝑛_𝑡𝑛−1𝐮_𝑡(𝑡)²₂𝑑𝑡. (97) We note 𝒩(2𝐮^𝑛_ℎ−𝐮_ℎ^𝑛−1,𝛿𝐄_ℎ^𝑛+1,𝛿𝐄_ℎ^𝑛+1)=0 which comes from (30). Then we obtain 𝐴_1,5+𝐴_1,6=−4𝜏𝒩2𝛿𝐮^𝑛_ℎ−𝛿𝐮_ℎ^𝑛−1,𝐆^𝑛+1+𝐄_ℎ^𝑛+1,𝛿𝐄_ℎ^𝑛+1−4𝜏𝒩2𝐮_ℎ^𝑛−1−𝐮_ℎ^𝑛−2,𝛿𝐆^𝑛+1,𝛿𝐄_ℎ^𝑛+1=4𝜏𝒩2𝛿𝐄^𝑛−𝛿𝐄^𝑛−1−2𝛿𝐮𝑡^𝑛+𝛿𝐮𝑡^𝑛−1,𝐆^𝑛+1+𝐄_ℎ^𝑛+1,𝛿𝐄_ℎ^𝑛+1+4𝜏𝒩2𝐄^𝑛−1−𝐄^𝑛−2−2𝐮𝑡^𝑛−1+𝐮𝑡^𝑛−2,𝛿𝐆^𝑛+1,𝛿𝐄_ℎ^𝑛+1∶=𝐵₁+𝐵₂. (98) To attack 𝐵₁, we first note Lemma 5 which is, for any 𝐰_ℎ∈𝕍_ℎ, ‖𝐰_ℎ‖_𝐋³_(Ω)≤𝐶ℎ^−𝑑/6‖𝐰_ℎ‖₀. If we apply 𝐄^𝑛+1+𝐆^𝑛+1₀+√𝜏+ℎ𝐄^𝑛+1+𝐆^𝑛+1₁≤𝐶(𝜏+ℎ) (99) which is the result of Lemma 13, then we can conclude, in light of (34), 𝐵₁≤𝐶𝜏2𝛿𝐄^𝑛−𝛿𝐄^𝑛−1_𝐋3(Ω)×𝐄_ℎ^𝑛+1+𝐆^𝑛+1₁∇𝛿𝐄_ℎ^𝑛+1₀+𝐶𝜏2𝛿𝐮(𝑡^𝑛)−𝛿𝐮(𝑡^𝑛−1)₂×𝐄_ℎ^𝑛+1+𝐆^𝑛+1₀∇𝛿𝐄_ℎ^𝑛+1₀≤𝜇𝜏6∇𝛿𝐄_ℎ^𝑛+12₀+𝐶𝜏𝜇𝜏1+ℎ×2𝛿𝐄^𝑛−𝛿𝐄^𝑛−12₀+𝐶𝜏²𝜇(𝜏+ℎ)^2𝑡𝑛_𝑡𝑛−2𝐮_𝑡(𝑡)²₂𝑑𝑡. (100) We now estimate 𝐵₂ using ‖2𝐄^𝑛−1−𝐄^𝑛−2‖_𝐋³_(Ω)≤𝐶ℎ^−𝑑/6‖2𝐄^𝑛−1−𝐄^𝑛−2‖₀≤𝑀: 𝐵₂≤𝐶𝜏2𝐄^𝑛−1−𝐄^𝑛−2_𝐋3(Ω)×𝛿𝐆^𝑛+1₁∇𝛿𝐄_ℎ^𝑛+1₀+𝐶𝜏2𝐮(𝑡^𝑛−1)−𝐮(𝑡^𝑛−2)₁×𝛿𝐆^𝑛+1₁∇𝛿𝐄_ℎ^𝑛+1₀≤𝐶𝜏𝛿𝐆^𝑛+12₁+𝜇𝜏6∇𝛿𝐄_ℎ^𝑛+12₀≤𝜇𝜏6∇𝛿𝐄_ℎ^𝑛+12₀+𝐶𝜏²ℎ²×^𝑡𝑛+1_𝑡𝑛−1𝐮_𝑡(𝑡)²₂+𝑝_𝑡(𝑡)²₁𝑑𝑡. (101) In light of Hölder inequality, (51) yields 𝐴₂=−23𝛿𝐆^𝑛+1−4𝛿𝐆^𝑛+𝛿𝐆^𝑛−1,𝛿𝐄_ℎ^𝑛+1≤𝜇𝜏6∇𝛿𝐄_ℎ^𝑛+12₀+𝐶ℎ^4𝑡𝑛+1_𝑡𝑛−2𝐮_𝑡(𝑡)²₂+𝑝_𝑡(𝑡)²₁𝑑𝑡. (102) Integral by parts leads to 𝐴₃=4𝜏𝛿𝛿𝑃_ℎ^𝑛+1,∇⋅𝛿𝐄_ℎ^𝑛+1≤𝜇𝜏6∇𝛿𝐄_ℎ^𝑛+12₀+𝐶𝜏⁴𝜇×^𝑡𝑛+1_𝑡𝑛−1𝐮_𝑡𝑡(𝑡)²₂+𝑝_𝑡𝑡(𝑡)²₁𝑑𝑡. (103) In order to tackle 𝐴₄, marking use of 𝛿𝐄_ℎ^𝑛+1=𝛿𝐄_ℎ^𝑛+1+∇𝛿𝛿𝜓_ℎ^𝑛+1, we readily get 𝐴₄=−4𝜏∇𝛿𝜀^𝑛_ℎ,𝛿𝐄_ℎ^𝑛+1+∇𝛿𝛿𝜓_ℎ^𝑛+1=−8𝜏²3∇𝛿𝜀^𝑛_ℎ,∇𝛿𝛿𝜀_ℎ^𝑛+1−𝛿𝛿𝑃_ℎ^𝑛+1≤−4𝜏²3∇𝛿𝜀_ℎ^𝑛+12₀−∇𝛿𝜀^𝑛_ℎ²₀−∇𝛿𝛿𝜀_ℎ^𝑛+12₀+𝐶𝜏³∇𝛿𝜀^𝑛_ℎ²₀+𝐶𝜏∇𝛿𝛿𝑃_ℎ^𝑛+12₀. (104) If we now apply inequality (𝑎+𝑏)²≤4𝑎²+(4/3)𝑏², then we can have 4𝜏²3∇𝛿𝛿𝜀_ℎ^𝑛+12₀=4𝜏²3∇𝛿𝛿𝑃_ℎ^𝑛+1+32𝜏∇𝛿𝛿𝜓_ℎ^𝑛+12₀≤𝐶𝜏²∇𝛿𝛿𝑃_ℎ^𝑛+12₀+4∇𝛿𝛿𝜓_ℎ^𝑛+12₀. (105) So we arrive at 𝐴₄≤−4𝜏²3∇𝛿𝜀_ℎ^𝑛+12₀−∇𝛿𝜀^𝑛_ℎ²₀+4∇𝛿𝛿𝜓_ℎ^𝑛+12₀+𝐶𝜏³∇𝛿𝜀^𝑛_ℎ²₀+𝐶𝜏^4𝑡𝑛+1_𝑡𝑛−1𝐮_𝑡𝑡(𝑡)²₂+𝑝_𝑡𝑡(𝑡)²₁𝑑𝑡. (106) Invoking (9), (37) and (60) lead to 𝐴₅=−4𝜇𝜏𝛿𝑞^𝑛_ℎ,𝛿𝛿𝑞_ℎ^𝑛+1=−2𝜇𝜏𝛿𝑞_ℎ^𝑛+12₀−𝛿𝑞^𝑛_ℎ²₀−𝛿𝛿𝑞_ℎ^𝑛+12₀≤−2𝜇𝜏𝛿𝑞_ℎ^𝑛+12₀−𝛿𝑞^𝑛_ℎ²₀+2𝜇𝜏∇𝛿𝐄_ℎ^𝑛+12₀. (107) Finally, the 𝐴₆ term becomes 𝐴₆≤𝐶𝜏𝛿𝐑^𝑛+1₀𝛿𝐄_ℎ^𝑛+1₀≤𝐶𝜏𝛿𝐄_ℎ^𝑛+12₀+𝐶𝜏^4𝑡𝑛+1_𝑡𝑛−1𝐮_𝑡𝑡𝑡(𝑡)²₀𝑑𝑡. (108) For the new term 𝐴₇, we argue as in Lemma 13 to arrive at 𝐴₇≤𝐶𝜅²𝜇⁴𝜏×2𝛿𝜗^𝑛_ℎ−𝛿𝜗_ℎ^𝑛−12₀+2𝛿𝜂^𝑛−𝛿𝜂^𝑛−12₀+𝜏𝛿𝐄_ℎ^𝑛+12₀+12∇𝛿𝛿𝜓_ℎ^𝑛+12₀+𝐶𝜏^4𝑡𝑛_𝑡⁺¹_𝑛−2𝜃_𝑡𝑡(𝑡)²₀𝑑𝑡. (109) We now insert the above estimates into (93) to obtain 𝛿𝛿𝐄_ℎ^𝑛+12₀+𝛿2𝛿𝐄_ℎ^𝑛+1−𝛿𝐄^𝑛_ℎ²₀+𝛿𝛿𝛿𝐄_ℎ^𝑛+12₀+𝜇𝜏∇𝛿𝐄_ℎ^𝑛+12₀+∇𝛿𝛿𝜓_ℎ^𝑛+12₀+2𝜇𝜏𝛿𝑞_ℎ^𝑛+12₀+4𝜏²3∇𝛿𝜀_ℎ^𝑛+12₀≤𝐶𝜏𝜇2𝛿𝐄^𝑛_ℎ−𝛿𝐄_ℎ^𝑛−12₀+2𝛿𝐆^𝑛−𝛿𝐆^𝑛−12₀+4𝜏²3∇𝛿𝜀^𝑛_ℎ²₀+𝐶𝜅²𝜇⁴𝜏×2𝛿𝜗^𝑛_ℎ−𝛿𝜗_ℎ^𝑛−12₀+2𝛿𝜂^𝑛−𝛿𝜂^𝑛−12₀+𝜏𝛿𝐄_ℎ^𝑛+12₀+2𝜇𝜏𝛿𝑞^𝑛_ℎ²₀+𝐶𝜏⁴×^𝑡𝑛+1_𝑡𝑛−1𝐮_𝑡𝑡𝑡(𝑡)²₀+𝐮_𝑡𝑡(𝑡)²₂+𝑝_𝑡𝑡(𝑡)²₁+𝜃_𝑡𝑡(𝑡)²₀𝑑𝑡+𝐶𝜏³∇𝛿𝜀^𝑛_ℎ²₀+𝐶𝜏²+ℎ²²×^𝑡𝑛+1_𝑡𝑛−2𝐮_𝑡𝑡(𝑡)²₀+𝐮_𝑡(𝑡)²₂+𝑝_𝑡(𝑡)²₁𝑑𝑡. (110)
To evaluate errors of 𝛿𝜗_ℎ^𝑛+1, we subtract two consecutive formulas (83) and choose by 𝜙_ℎ=4𝜏𝛿𝜗_ℎ^𝑛+1=4𝜏(𝛿𝜗^𝑛+1−𝛿𝜂^𝑛+1) as follows: 𝛿𝛿𝜗_ℎ^𝑛+12₀+𝛿2𝛿𝜗_ℎ^𝑛+12₀+𝛿𝛿𝛿𝜗_ℎ^𝑛+12₀+4𝜆𝜇𝜏∇𝛿𝜗_ℎ^𝑛+12₀=𝐴₈+𝐴₉, (111) where 𝐴₈∶=4𝜏𝒩𝐮𝑡^𝑛,𝜃𝑡^𝑛,𝛿𝜗_ℎ^𝑛+1−4𝜏𝒩𝐮^𝑛_ℎ,𝜃^𝑛_ℎ,𝛿𝜗_ℎ^𝑛+1−4𝜏𝒩𝐮𝑡^𝑛+1,𝜃𝑡^𝑛+1,𝛿𝜗_ℎ^𝑛+1+4𝜏𝒩𝐮_ℎ^𝑛+1,𝜃_ℎ^𝑛+1,𝛿𝜗_ℎ^𝑛+1,𝐴₉∶=−23𝛿𝜂^𝑛+1−4𝛿𝜂^𝑛+𝛿𝜂^𝑛−1,𝛿𝜗_ℎ^𝑛+1+4𝜏𝛿𝑄^𝑛+1,𝛿𝜗_ℎ^𝑛+1. (112) We treat the 𝐴₈ term first by rewriting it as follows: 𝐴₈=−4𝜏𝒩𝐄^𝑛+1,𝜃𝑡^𝑛+1,𝛿𝜗_ℎ^𝑛+1−4𝜏𝒩𝐮_ℎ^𝑛+1,𝜗^𝑛+1,𝛿𝜗_ℎ^𝑛+1+4𝜏𝒩𝐄^𝑛,𝜃𝑡^𝑛,𝛿𝜗_ℎ^𝑛+1+4𝜏𝒩𝐮^𝑛_ℎ,𝜗^𝑛,𝛿𝜗_ℎ^𝑛+1∶=⁴_𝑖=1𝐴_8,𝑖. (113) In estimating convection terms, we will use Lemma 10 frequently without notice. We recall ‖𝜃(𝑡)‖₂≤𝐶 and ‖𝐄^𝑛‖₀≤𝐶(𝜏+ℎ) which is the result in Lemma 13 to obtain 𝐴_8,1+𝐴_8,3=−4𝜏𝒩𝛿𝐄^𝑛+1,𝜃𝑡^𝑛+1,𝛿𝜗_ℎ^𝑛+1−4𝜏𝒩𝐄^𝑛,𝛿𝜃𝑡^𝑛+1,𝛿𝜗_ℎ^𝑛+1≤𝐶𝜏𝛿𝐄^𝑛+1₀𝜃(𝑡^𝑛+1)₂+𝐄^𝑛₀𝛿𝜃(𝑡^𝑛+1)₂𝛿𝜗_ℎ^𝑛+1₁≤𝜆𝜇𝜏∇𝛿𝜗_ℎ^𝑛+12₀+𝐶𝜏𝜆𝜇𝛿𝐄^𝑛+12₀+𝐶𝜏²(𝜏+ℎ)²×𝜆𝜇_{𝑡𝑛+1𝑡𝑛}𝜃_𝑡(𝑡)²₂𝑑𝑡. (114) In conjunction with 𝒩(𝐮_ℎ^𝑛+1,𝛿𝜗_ℎ^𝑛+1,𝛿𝜗_ℎ^𝑛+1)=0 which comes from (30), we rewrite 𝐴_8,2+𝐴_8,4 as 𝐴_8,2+𝐴_8,4=4𝜏𝒩𝐄_ℎ^𝑛+1−𝐔_ℎ^𝑛+1,𝛿𝜂^𝑛+1,𝛿𝜗_ℎ^𝑛+1+4𝜏𝒩𝛿𝐄_ℎ^𝑛+1−𝛿𝐔_ℎ^𝑛+1,𝜗^𝑛,𝛿𝜗_ℎ^𝑛+1∶=𝐵₃+𝐵₄. (115) Because we can readily get ‖𝐄_ℎ^𝑛+1−𝐔_ℎ^𝑛+1‖_𝐋³_(Ω)≤𝑀 via Lemmas 5 and 13, we conclude, in conjunction with (60), 𝐵₃≤𝐶𝜏𝐄_ℎ^𝑛+1−𝐔_ℎ^𝑛+1_𝐋3(Ω)×𝛿𝜂^𝑛+1₁𝛿𝜗_ℎ^𝑛+1₁≤𝜆𝜇𝜏∇𝛿𝜗_ℎ^𝑛+12₀+𝐶𝜏²ℎ²×𝜆𝜇_{𝑡𝑛+1𝑡𝑛}𝜃_𝑡(𝑡)²₂𝑑𝑡. (116) Before we attack 𝐵₄, we note that the assumption 𝜏≤𝐶ℎ is required to apply ‖𝛿𝐄^𝑛_ℎ‖_𝐋³_(Ω)≤𝐶𝜏^−1/2‖𝛿𝐄^𝑛_ℎ‖₀. We now imply ‖𝜗^𝑛‖²₀+𝜏‖𝜗^𝑛‖²₁≤𝐶(𝜏²+ℎ²) and (51) to get 𝐵₄=4𝜏𝒩𝛿𝐆^𝑛+1−𝛿𝐮𝑡^𝑛+1+𝛿𝐄_ℎ^𝑛+1,𝜗^𝑛,𝛿𝜗_ℎ^𝑛+1≤𝐶𝜏𝛿𝐆^𝑛+1₁+𝛿𝐄_ℎ^𝑛+1_𝐋3(Ω)×𝜗^𝑛₁+𝛿𝐮(𝑡^𝑛)₂𝜗^𝑛₀𝛿𝜗_ℎ^𝑛+1₁≤𝐶𝜏𝜆𝜇𝛿𝐄_ℎ^𝑛+12₀+𝜆𝜇𝜏∇𝛿𝜗_ℎ^𝑛+12₀+𝐶𝜏²𝜏²+ℎ²×𝜆𝜇_{𝑡𝑛+1𝑡𝑛−1}𝐮_𝑡(𝑡)²₂+𝑝_𝑡(𝑡)²₁𝑑𝑡. (117) In light of Hölder inequality, 𝐴₉ becomes 𝐴₉≤C𝜏𝛿𝜗_ℎ^𝑛+12₀+𝐶𝜏⁴+ℎ⁴×^𝑡𝑛+1_𝑡𝑛−2𝜃_𝑡(𝑡)²₂+𝜃_𝑡𝑡𝑡²₀𝑑𝑡. (118) Inserting the above estimates into (111) yields 𝛿𝛿𝜗_ℎ^𝑛+12₀+𝛿2𝛿𝜗_ℎ^𝑛+1−𝛿𝜗^𝑛_ℎ²₀+𝛿𝛿𝛿𝜗_ℎ^𝑛+12₀+𝜆𝜇𝜏∇𝛿𝜗_ℎ^𝑛+12₀≤𝐶𝜏𝛿𝜗_ℎ^𝑛+12₀+𝐶𝜏𝜆𝜇𝛿𝐄_ℎ^𝑛+1₀+𝛿𝐄^𝑛_ℎ²₀+𝛿𝐆^𝑛+1₀+𝐶𝜏²+ℎ²²×^𝑡𝑛+1_𝑡𝑛−2𝐮_𝑡(𝑡)²₂+𝑝_𝑡(𝑡)²₁+𝜃_𝑡(𝑡)²₂+𝜃_𝑡𝑡𝑡²₀𝑑𝑡. (119) Adding 2 equations (110) and (119) and then summing up 𝑛 from 2 to 𝑁 lead up to (92) and complete the proof.

We now estimate optimal accuracy for velocity and temperature.

Lemma 17 (full rate of convergence for velocity and temperature). One denotes that (𝐯^𝑛+1,𝑟^𝑛+1) and (𝐯_ℎ^𝑛+1,𝑟_ℎ^𝑛+1) are the solutions of (16) and (24) with 𝐳=𝐄_ℎ^𝑛+1, respectively. And let 𝜔^𝑛+1 and 𝜔_ℎ^𝑛+1 be solutions of (17) and (28) with 𝜉=𝜗_ℎ^𝑛+1. Let the exact solution of (1) be smooth enough and 𝜏≤𝐶ℎ. If Assumptions 1 and 3–5 hold, then one has ∇𝐯_ℎ^𝑁+12₀+∇2𝐯_ℎ^𝑁+1−𝐯^𝑁_ℎ²₀+∇𝜔_ℎ^𝑁+12₀+∇2𝜔_ℎ^𝑁+1−𝜔^𝑁_ℎ²₀+^𝑁_𝑛=1∇𝛿𝛿𝐯_ℎ^𝑛+12₀+∇𝛿𝛿𝜔_ℎ^𝑛+12₀+2𝜇𝜏^𝑁_𝑛=1𝐄_ℎ^𝑛+12₀+𝜆𝜗_ℎ^𝑛+12₀≤𝐶𝜏⁴+ℎ⁴. (120)

Proof. We choose 𝐰_ℎ=4𝜏𝐯_ℎ^𝑛+1∈𝕍_ℎ in (66) to obtain 𝛿∇𝐯_ℎ^𝑛+12₀+𝛿∇2𝐯_ℎ^𝑛+1−𝐯^𝑛_ℎ²₀+∇𝛿𝛿𝐯_ℎ^𝑛+12₀+4𝜇𝜏𝐄_ℎ^𝑛+12₀=⁵_𝑖=1𝐴_𝑖, (121) where 𝐴₁∶=4𝜏𝒩2𝐮^𝑛_ℎ−𝐮_ℎ^𝑛−1,𝐮_ℎ^𝑛+1,𝐯_ℎ^𝑛+1−4𝜏𝒩𝐮𝑡^𝑛+1,𝐮𝑡^𝑛+1,𝐯_ℎ^𝑛+1,𝐴₂∶=3𝐆^𝑛+1−4𝐆^𝑛+𝐆^𝑛−1,𝐯_ℎ^𝑛+1,𝐴₃∶=4𝜇𝜏𝐄_ℎ^𝑛+1,∇𝑟_ℎ^𝑛+1,𝐴₄∶=4𝜏𝐑^𝑛+1,𝐯_ℎ^𝑛+1,𝐴₅∶=4𝜅𝜇²𝜏𝐠𝜃𝑡^𝑛+1−2𝜃^𝑛_ℎ+𝜃_ℎ^𝑛−1,𝐯_ℎ^𝑛+1. (122) We now estimate all the terms from 𝐴₁ to 𝐴₅, respectively. The convection term 𝐴₁ can be rewritten as follows: 𝐴₁=−4𝜏𝒩𝛿𝛿𝐮𝑡^𝑛+1,𝐮𝑡^𝑛+1,𝐯_ℎ^𝑛+1−4𝜏𝒩2𝐮𝑡^𝑛−𝐮𝑡^𝑛−1,𝐄^𝑛+1,𝐯_ℎ^𝑛+1+4𝜏𝒩2𝐄^𝑛−𝐄^𝑛−1,𝐄^𝑛+1,𝐯_ℎ^𝑛+1−4𝜏𝒩2𝐄^𝑛−𝐄^𝑛−1,𝐮𝑡^𝑛+1,𝐯_ℎ^𝑛+1, (123) and we denote by 𝐴_1,𝑖, for 𝑖=1,2,…,4, the four terms in the right-hand side. To estimate convection terms, we will use frequently Lemma 10 without notice. Using ‖𝐮(𝑡^𝑛+1)‖₂≤𝑀, we can readily get 𝐴_1,1≤𝐶𝜏𝛿𝛿𝐮(𝑡^𝑛+1)₀×𝐮(𝑡^𝑛+1)₂𝐯_ℎ^𝑛+1₁≤𝐶𝜏∇𝐯_ℎ^𝑛+12₀+𝐶𝜇𝜏⁴×^𝑡𝑛+1_𝑡𝑛−1𝐮_𝑡𝑡(𝑡)²₀𝐴𝑑𝑡,_1,4≤𝐶𝜏𝐄^𝑛₀+𝛿𝐄^𝑛₀×𝐮(𝑡^𝑛+1)₂𝐯_ℎ^𝑛+1₁≤𝜇𝜏2𝐄^𝑛_ℎ²₀+𝐆^𝑛2₀+𝛿𝐄^𝑛_ℎ²₀+𝛿𝐆^𝑛2₀+𝐶𝜏𝜇∇𝐯_ℎ^𝑛+12₀. (124)(125) Because ∇⋅(2𝐮(𝑡^𝑛)−𝐮(𝑡^𝑛−1))=0 and 2𝐮(𝑡^𝑛)−𝐮(𝑡^𝑛−1)=𝟎 on boundary, we can use (32) and so we get 𝐴_1,2≤𝐶𝜏2𝐮(𝑡^𝑛)−𝐮(𝑡^𝑛−1)₂×𝐄^𝑛+1₀𝐯_ℎ^𝑛+1₁≤𝜇𝜏2𝐄_ℎ^𝑛+12₀+𝐆^𝑛+12₀+∇𝛿𝜓_ℎ^𝑛+12₀+𝐶𝜏𝜇∇𝐯_ℎ^𝑛+12₀. (126) In light of 𝜀_ℎ^𝑛+1=𝑃_ℎ^𝑛+1+3𝜓_ℎ^𝑛+1/2𝜏, we can obtain ∇𝛿𝜓_ℎ^𝑛+12₀=4𝜏²9∇𝛿𝜀_ℎ^𝑛+1−𝛿𝑃_ℎ^𝑛+12₀≤𝐶𝜏²∇𝛿𝜀_ℎ^𝑛+12₀+𝐶𝜏^3𝑡𝑛+1_𝑡𝑛𝐮_𝑡(𝑡)²₂+𝑝_𝑡(𝑡)²₁𝑑𝑡, (127) and so we can conclude 𝐴_1,2≤𝜇𝜏2𝐄_ℎ^𝑛+12₀+𝐆^𝑛+12₀+𝐶𝜇𝜏³∇𝛿𝜀_ℎ^𝑛+12₀+𝐶𝜇𝜏⁴×^𝑡𝑛+1_𝑡𝑛𝐮_𝑡(𝑡)²₂+𝑝_𝑡(𝑡)²₁+𝑑𝑡𝐶𝜏𝜇∇𝐯_ℎ^𝑛+12₀. (128) If we apply ‖2𝐄^𝑛−𝐄^𝑛−1‖_𝐋³_(Ω)≤𝐶ℎ^−𝑑/6‖2𝐄^𝑛−𝐄^𝑛−1‖₀≤𝐶ℎ^−𝑑/6(𝜏+ℎ) which can be derived from Lemmas 5 and 13, then we can get, by the help of Lemma 10 and Assumption 1, 𝐴_1,3≤𝐶𝜏2𝐄^𝑛−𝐄^𝑛−1₀×𝐄^𝑛+1₁𝐯^𝑛+1₂+𝐶𝜏2𝐄^𝑛−𝐄^𝑛−1_𝐋3(Ω)×𝐄^𝑛+1₁𝐯^𝑛+1−𝐯_ℎ^𝑛+1₁≤𝐶𝜏(𝜏+ℎ)𝐄^𝑛+1₁𝐯^𝑛+1₂+𝐶𝜏ℎ𝐄^𝑛+1₁𝐯^𝑛+1₂≤𝐶𝜇𝜏(𝜏+ℎ)²×∇𝐄_ℎ^𝑛+12₀+∇𝐆^𝑛+12₀+𝜇𝜏2𝐄_ℎ^𝑛+12₀. (129) In conjunction with (51), we can have 𝐴₂≤𝐶𝜏𝛿𝐆^𝑛+12₀+𝛿𝐆^𝑛2₀+𝐶𝜏𝐯_ℎ^𝑛+12₀≤𝐶𝜏𝐯_ℎ^𝑛+12₀+𝐶ℎ⁴×^𝑡𝑛+1_𝑡𝑛−1𝐮_𝑡(𝑡)²₂+𝑝_𝑡(𝑡)²₁𝑑𝑡. (130) The definition 𝐄_ℎ^𝑛+1=𝐄_ℎ^𝑛+1+∇𝛿𝜓_ℎ^𝑛+1 and Assumption 1 give us 𝐴₃=4𝜇𝜏∇𝛿𝜓_ℎ^𝑛+1,∇𝑟_ℎ^𝑛+1≤𝜇𝜏𝐄_ℎ^𝑛+12₀+𝐶𝜏𝜇∇𝛿𝜓_ℎ^𝑛+12₀. (131) If we apply (127) again, then we arrive at 𝐴₃≤𝜇𝜏𝐄_ℎ^𝑛+12₀+𝐶𝜏³𝜇∇𝛿𝜀_ℎ^𝑛+12₀+𝐶𝜏⁴𝜇^𝑡𝑛+1_𝑡^𝑛∇𝑝_𝑡(𝑡)²₀𝑑𝑡. (132) On the other hand, the truncation error term becomes 𝐴₄=4𝜏𝐑^𝑛+1,𝐯_ℎ^𝑛+1≤𝐶𝜏∇𝐯_ℎ^𝑛+12₀+𝐶𝜏⁴×^𝑡𝑛+1_𝑡𝑛−1𝐮_𝑡𝑡𝑡(𝑠)²₀𝑑𝑡. (133) For the new term 𝐴₅, we observe ‖𝛿𝛿𝜃(𝑡^𝑛+1)‖₀≤𝐶𝜏³∫^𝑡𝑛+1_𝑡𝑛−1‖𝜃_𝑡(𝑡)‖²₀𝑑𝑡; whence 𝐴₅=4𝜅𝜇²𝜏×𝐠𝛿𝛿𝜃𝑡^𝑛+1+2𝜗^𝑛−𝜗^𝑛−1,𝐯_ℎ^𝑛+1≤𝐶𝜏∇𝐯_ℎ^𝑛+12₀+𝐶𝜅²𝜇⁴𝜏×2𝜗^𝑛−𝜗^𝑛−12₀+𝐶𝜅²𝜇⁴𝜏⁴×^𝑡𝑛+1_𝑡𝑛−1𝜃_𝑡(𝑡)²₀𝑑𝑡. (134) Invoking 𝐯⁰=𝟎 and inserting the above estimates from 𝐴₁ and 𝐴₅ into (121) lead us to 𝛿∇𝐯_ℎ^𝑛+12₀+𝛿∇2𝐯_ℎ^𝑛+1−𝐯^𝑛_ℎ²₀+∇𝛿𝛿𝐯_ℎ^𝑛+12₀+2𝜇𝜏𝐄_ℎ^𝑛+12₀≤𝜇𝜏2𝐄^𝑛_ℎ²₀+𝐆^𝑛+12₀+𝐆^𝑛2₀+𝛿𝐄^𝑛_ℎ²₀+𝛿𝐆^𝑛2₀+𝐶𝜅²𝜇⁴𝜏2𝜗^𝑛−𝜗^𝑛−12₀+𝐶𝜇𝜏(𝜏+ℎ)²×∇𝐄_ℎ^𝑛+12₀+∇𝐆^𝑛+12₀+𝐶𝜏³𝜇∇𝛿𝜀_ℎ^𝑛+12₀+𝐶𝜇𝜏⁴+ℎ⁴×^𝑡𝑛+1_𝑡𝑛−1𝐮_𝑡(𝑡)²₂+𝑝_𝑡(𝑡)²₁𝑑𝑡+𝐶𝜏∇𝐯_ℎ^𝑛+12₀+𝐶𝜏⁴×^𝑡𝑛+1_𝑡𝑛−1𝐮_𝑡𝑡(𝑡)²₀+𝐮_𝑡𝑡𝑡(𝑡)²₀+𝜃_𝑡(𝑡)²₀𝑑𝑡. (135) On the other hand, we choose 𝜙_ℎ=4𝜏𝜔_ℎ^𝑛+1 in (83) to obtain 𝛿∇𝜔_ℎ^𝑛+12₀+𝛿∇2𝜔_ℎ^𝑛+1−𝜔^𝑛_ℎ²₀+∇𝛿𝛿𝜔_ℎ^𝑛+12₀+4𝜆𝜇𝜏𝜗_ℎ^𝑛+12₀=𝐴₆+𝐴₇, (136) where 𝐴₆∶=4𝜏𝒩𝐮_ℎ^𝑛+1,𝜃_ℎ^𝑛+1,𝜔_ℎ^𝑛+1−4𝜏𝒩𝐮𝑡^𝑛+1,𝜃𝑡^𝑛+1,𝜔_ℎ^𝑛+1,𝐴₇∶=−23𝜂^𝑛+1−4𝜂^𝑛+𝜂^𝑛−1,𝜔_ℎ^𝑛+1+4𝜏𝑄^𝑛+1,𝜔_ℎ^𝑛+1. (137) In order to estimate 𝐴₆, we note first ‖𝐄^𝑛+1‖_𝐋³_(Ω)≤𝐶ℎ^−1/2‖𝐄^𝑛+1‖₀ and ‖𝜗^𝑛+1‖²₁≤𝐶((𝜏²+ℎ²)/𝜏) which come from Lemmas 5 and 13, respectively. Then Lemma 10 and the assumption 𝜏≤𝐶ℎ yield 𝐴₆=−4𝜏𝒩𝐄^𝑛+1,𝜃𝑡^𝑛+1,𝜔_ℎ^𝑛+1+4𝜏𝒩𝐄^𝑛+1−𝐮𝑡^𝑛+1,𝜗^𝑛+1,𝜔_ℎ^𝑛+1≤𝐶𝜏𝐄^𝑛+1₀𝜃(𝑡^𝑛+1)₂𝜔_ℎ^𝑛+1₁+𝐶𝜏𝐄^𝑛+1_𝐋3(Ω)𝜗^𝑛+1₁𝜔_ℎ^𝑛+1₁+𝐶𝜏𝐮(𝑡^𝑛+1)₂𝜗^𝑛+1₀𝜔_ℎ^𝑛+1₁≤𝐶𝜏∇𝜔_ℎ^𝑛+12₀+𝜇𝜏𝐄^𝑛+12₀+𝜆𝜇𝜏𝜗_ℎ^𝑛+12₀+𝜂^𝑛+12₀. (138) Since ‖𝛿𝜂^𝑛+1‖²₀≤𝐶ℎ⁴𝜏∫_𝑛^𝑛+1‖𝜃_𝑡(𝑡)‖²₂𝑑𝑡, we can readily get 𝐴₇≤𝐶𝜏∇𝜔_ℎ^𝑛+12₀+𝐶𝜏⁴+ℎ⁴×^𝑛+1_𝑛−1𝜃_𝑡(𝑡)²₂+𝜃_𝑡𝑡𝑡(𝑡)²₂𝑑𝑡. (139) Inserting the above estimates into (136) yields 𝛿∇𝜔_ℎ^𝑛+12₀+𝛿∇2𝜔_ℎ^𝑛+1−𝜔^𝑛_ℎ²₀+∇𝛿𝛿𝜔_ℎ^𝑛+12₀+3𝜆𝜇𝜏𝜗_ℎ^𝑛+12₀≤𝐶𝜏∇𝜔_ℎ^𝑛+12₀+𝜆𝜇𝜏𝜂^𝑛+12₀+𝜇𝜏𝐄^𝑛+12₀+𝐶𝜏⁴+ℎ⁴×^𝑛+1_𝑛−1𝜃_𝑡(𝑡)²₂+𝜃_𝑡𝑡𝑡(𝑡)²₂𝑑𝑡. (140)
In conjunction with the discrete Gronwall inequality, adding (135) and (140) and summing over 𝑛 from 1 to 𝑁 lead to (120).

We now estimate the pressure error in 𝐿²(0,𝑇;𝐿²(Ω)). This hinges on the error estimates for the time derivative of velocity and temperature of Lemma 15.

Lemma 18 (pressure error estimate). Let the exact solution of (1) be smooth enough and 𝜏≤𝐶ℎ. If Assumptions 1–5 hold, then one has 𝜏^𝑁_𝑛=1𝑒_ℎ^𝑛+12₀≤𝐶𝜏²+ℎ². (141)

Proof. We first recall again inf-sup condition in Assumption 2. Consequently, it suffices to estimate ⟨𝑒^𝑛+1,∇⋅𝐰⟩ in terms of ‖∇𝐰‖₀. In conjunction with (10), we can rewrite (66) as 𝑒_ℎ^𝑛+1,∇⋅𝐰_ℎ=12𝜏3𝐄^𝑛+1−4𝐄^𝑛+𝐄^𝑛−1,𝐰_ℎ+𝜇∇𝐄_ℎ^𝑛+1,∇𝐰_ℎ+𝒩𝛿𝛿𝐮𝑡^𝑛+1,𝐮𝑡^𝑛+1,𝐰_ℎ+𝒩2𝐮𝑡^𝑛+1−𝐮𝑡^𝑛,𝐄^𝑛+1,𝐰_ℎ+𝒩2𝐄^𝑛−𝐄^𝑛−1,𝐮_ℎ^𝑛+1,𝐰_ℎ−𝜇∇𝛿𝑞_ℎ^𝑛+1,𝐰_ℎ−𝐑^𝑛+1,𝐰_ℎ+𝜅𝜇²𝐠𝜃𝑡^𝑛+1−2𝜃^𝑛_ℎ+𝜃_ℎ^𝑛−1,𝐰_ℎ∶=⁸_𝑖=1𝐴_𝑖. (142) We now proceed to estimate each term from 𝐴₁ to 𝐴₇ separately. We readily obtain 𝐴₁≤𝐶𝜏𝛿𝐄^𝑛+1₀+𝛿𝐄^𝑛₀𝐰_ℎ0≤𝐶𝜏𝛿𝐄^𝑛+1₀+𝛿𝐄^𝑛₀∇𝐰_ℎ0,𝐴₂≤𝐶∇𝐄_ℎ^𝑛+1₀∇𝐰_ℎ0. (143) Terms 𝐴₃ and 𝐴₄ can be dealt with thanks to the aid of Lemma 10 and ‖𝐮(𝑡^𝑛+1)‖₂≤𝑀 as follows: 𝐴₃≤𝐶𝛿𝛿𝐮𝑡^𝑛+1₀×𝐮(𝑡^𝑛+1)₂𝐰_ℎ1≤𝐶𝛿𝛿𝐮𝑡^𝑛+1₁∇𝐰_ℎ1,𝐴₄≤𝐶2𝐮𝑡^𝑛+1−𝐮𝑡^𝑛₂×𝐄^𝑛+1₀𝐰_ℎ1≤𝐶𝐄^𝑛+1₀∇𝐰_ℎ0. (144) In light of ‖𝐮_ℎ^𝑛+1‖₁=‖𝐄^𝑛+1−𝐮(𝑡^𝑛+1)‖₁≤𝐶 from Lemma 13, we can have 𝐴₅≤𝐶2𝐄^𝑛−𝐄^𝑛−1_𝐋3(Ω)×𝐮_ℎ^𝑛+1₁𝐰_ℎ1≤𝐶√ℎ𝐄^𝑛+1₀+𝐄^𝑛₀×∇𝐰_ℎ0. (145) Integrate by parts and Hölder inequality yield 𝐴₆≤𝐶𝛿𝑞_ℎ^𝑛+1₀∇𝐰_ℎ0. (146) On the other hand, we have 𝐴₇≤𝐑^𝑛+1₋₁‖∇𝐰‖₀. (147) The new term 𝐴₈ can be bound by the Hölder inequality as 𝐴₈=𝜅𝜇²×𝐠𝛿𝛿𝜃𝑡^𝑛+1−2𝜗^𝑛+𝜗^𝑛−1,𝐰_ℎ^𝑛+1≤𝐶𝜅𝜇²𝛿𝛿𝜃(𝑡^𝑛+1)₀+2𝜗^𝑛_ℎ−𝜗_ℎ^𝑛−1₀+2𝜂^𝑛−𝜂^𝑛−1₀‖∇𝐰‖₀. (148) Inserting the estimates from 𝐴₁ to 𝐴₇ back into (142) and employing the discrete inf-sup condition in Assumption 2, we obtain 𝐶𝑒^𝑛+1₀≤1𝜏𝛿𝐄^𝑛+1₀+𝛿𝐄^𝑛₀+1√ℎ𝐄^𝑛+1₀+𝐄^𝑛₀+𝛿𝑞_ℎ^𝑛+1₀+∇𝐄_ℎ^𝑛+1₀+𝛿𝛿𝐮(𝑡^𝑛+1)₁+𝐄^𝑛+1₀+𝑅^𝑛+1₋₁+𝛿𝛿𝜃(𝑡^𝑛+1)₀+2𝜗^𝑛_ℎ−𝜗_ℎ^𝑛−1₀+2𝜂^𝑛−𝜂^𝑛−1₀. (149) If we now square it, multiply it by 𝜏, and sum over 𝑛 from 1 to 𝑁, then Lemmas 13, 15, and 17 derive (141).

5. Numerical Experiments

We finally document 3 computational performances of SGUM. The first is to check accuracy, and then the next 2 examples are physically relevant numerical simulations, the Bénard convection problem and the thermal driven cavity flow. We perform the last 2 examples under the same set within [2], but we conclude with different numerical simulation for the second test, the Bénard convection problem, from that of [2]. We impose Taylor-Hood (𝒫2−𝒫1) in all 3 experiments.

Example 19 (mesh analysis). In this first experiment, we choose square domain [0,1]×[0,1] and impose forcing term the exact solution to become 𝑢=𝜋𝑡²−𝑡+1𝑥²(1−𝑥)²,𝑡sin2𝜋𝑦𝑣=−2²−𝑡+1𝑥(2𝑥−1)(𝑥−1)sin²,𝑡𝜋𝑦𝑝=−²−𝑡−1cos(𝜋𝑥)𝑦²,(𝑡)()𝑦.+1𝜃=cossin𝜋𝑥1−𝑦 (150) Table 1 is error decay with 𝜇=1 and 𝜏=ℎ. We conclude that the numerical accuracy of SGUM is optimal and consists with the result of Theorem 4.

Table 1 

Error decay for Algorithm 1 with 𝜏=ℎ and 𝜇=1.

Example 20 (Bénard convection). In order to explore the applicability of the SGUM, we consider the Bénard convection on the domain Ω=[0,5]×[0,1] with forcing 𝐟=0 and 𝑏=0. Figure 1 displays the initial and boundary conditions for velocity 𝐮 and temperature 𝜃, as already studied in [2]. Figures 2–4 are simulations at 𝑡=1 with the nondimensional parameters 𝜅=10⁴,𝜆=1, and  𝜇=1. Figure 2 is the result for the case 𝜏=10⁻²,ℎ=2⁻⁴ which is the same condition in [2], and so it displays similar behavior within [2] including 6 circulations in the velocity stream line. However, Figures 3 and 4, the higher resolution simulations with 𝜏=10⁻³,ℎ=2⁻⁴, and 𝜏=10⁻⁴,ℎ=2⁻⁶, display 8 circulations in the stream line. So we conclude that the high resolution result is correct simulation and thus Figure 2 and the result in [2] are not eventual simulation.

Figure 1 

Example 20: initial and boundary values of Bénard convection problem.

Figure 2 

Example 20: streamlines of velocity and isolines of temperature and pressure, at time 𝑡=1.0. The nondimensional parameters are 𝜅=104,  𝜆=1,  𝜇=1, and the discretization parameters are 𝜏=10−2,  ℎ=2−4.

Figure 3 

Example 20: streamlines of velocity and isolines of temperature and pressure, at time 𝑡=1.0. The nondimensional parameters are 𝜅=104,𝜆=1, and 𝜇=1, and the discretization parameters are 𝜏=10−3,ℎ=2−4.

Figure 4 

Example 20: streamlines of velocity and isolines of temperature and pressure, at time 𝑡=1.0. The nondimensional parameters are 𝜅=104,𝜆=1, and 𝜇=1, and the discretization parameters are 𝜏=10−4,ℎ=2−6.

Example 21 (thermal driven cavity flow). We consider the thermal driven cavity flow in an enclosed square Ω=[0,1]², as already studied in several papers [2, 6, 7]. The experiment is carried out with the same setting as in the work of Gresho et al. [6], which is shown in Figure 5. Figure 6 displays the evolution from rest (𝑡=0) to steady state (𝑡=0.2).

Figure 5 

Example 21: initial and boundary values for thermal driven cavity flow.

(a)

(b)

(c)

(d)

Figure 6 

Example 21: time sequence 𝑡=0.003, 0.01, 0.02, 0.04, and 0.1 for the driven cavity. The first two columns are the streamlines and vector fields for velocity, and the third and fourth ones are the contour lines for pressure and temperature, respectively. The nondimensional parameters are 𝜅=105,𝜆=1, and 𝜇=1, and the discretization parameters are 𝜏=10−4,ℎ=2−5. Note that 𝑢max stands for ‖𝐮‖∞.