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

Authors: Jae-Hong Pyo
Journal: Journal of Applied Mathematics
 (2013)

---

## 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.


---

## Body

## 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 ℝd with d=2 or 3,
(1)ut+(u·∇)u+∇p-μΔu+κμ2gθ=f,inΩ,∇·u=0,inΩ,θt+u·∇θ-λμΔθ=b,inΩ,
with initial conditions u(x,0)=u0,  θ(x,0)=θ0 in Ω, vanishing Dirichlet boundary conditions u=0,  θ=0 on ∂Ω, and pressure mean-value ∫Ω‍p=0. The forcing functions f and b are given, and g is the vector of gravitational acceleration. The nondimensional numbers μ=Re-1 and λ=Pr-1 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 κμ2gθ 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 H1(Ω). 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 zn+1,
(2)δzn+1=zn+1-zn,δδzn+1=δ(δzn+1)=zn+1-2zn+zn-1,….In order to introduce the finite element discretization, we need further notations. LetHs(Ω) be the Sobolev space with s derivatives in L2(Ω), set L2(Ω)=(L2(Ω))d and Hs(Ω)=(Hs(Ω))d, where d=2 or 3, and denote by L02(Ω) the subspace of L2(Ω) of functions with vanishing mean value. We indicate with ∥·∥s the norm in Hs(Ω) and with 〈·,·〉 the inner product in L2(Ω). Let 𝔗={K} be a shape-regular quasiuniform partition of Ω of mesh size h into closed elements K [8–10]. The vector and scalar finite element spaces are
(3)𝕎h:={wh∈L2(Ω):wh|K∈𝒫(K)∀K∈𝔗},𝕍h:=𝕎h∩H01(Ω),𝕋h:={θh∈H01(Ω)∩C0(Ω):θh|K∈𝒬(K)∀K∈𝔗},ℙh:={qh∈L02(Ω)∩C0(Ω):qh|K∈ℛ(K)∀K∈𝔗},
where 𝒫(K), 𝒬(K), and ℛ(K) are spaces of polynomials with degree bounded uniformly with respect to K∈𝔗 [9, 10]. We stress that the space ℙh is composed of continuous functions to ensure the crucial equality
(4)〈∇·wh,ϕh〉=-〈wh,∇ϕh〉,∀wh∈𝕍h,∀ϕh∈ℙh.
Using the following discrete counterpart of the form 𝒩(u,v,w):=〈(u·∇)v,w〉(5)𝒩(uh,vh,wh):=12〈(uh·∇)vh,wh〉-12〈(uh·∇)wh,vh〉,
we now introduce the fully discrete SGUM to solve the evolution Boussinesq equations (1).Algorithm 1 (the fully discrete stabilized Gauge-Uzawa method).
Computeθh1, uh1, and ph1 via any method satisfying Assumption 4 below, and set ψh1=(-2τ/3)ph1 and qh1=0. Repeat for 1≤n≤N=[T/τ-1].Step 1.
Setuh*=2uhn-uhn-1 and θh*=2θhn-θhn-1, and then find u^hn+1∈𝕍h as the solution of
(6)12τ〈3u^hn+1-4uhn+uhn-1,wh〉+〈∇phn,wh〉+𝒩(uh*,u^hn+1,wh)+μ〈∇u^hn+1,∇wh〉+κμ2〈gθh*,wh〉=〈f(tn+1),wh〉,∀wh∈𝕍h.Step 2.
Findψhn+1∈ℙh as the solution of
(7)〈∇ψhn+1,∇ϕh〉=〈∇ψhn,∇ϕh〉+〈∇·u^hn+1,ϕh〉,〈∇ψhn,∇ϕh〉88.88.+∀ϕh∈ℙh.Step 3.
Updateuhn+1 and qhn+1∈ℙh according to
(8)uhn+1=u^hn+1+∇(ψhn+1-ψhn),(9)〈qhn+1,ϕh〉=〈qhn,ϕh〉-〈∇·u^hn+1,ϕh〉,∀ϕh∈ℙh.Step 4.
Update pressurephn+1 by
(10)phn+1=-3ψhn+12τ+μqhn+1.Step 5.
Findθhn+1∈𝕋h as the solution of
(11)12τ〈3θhn+1-4θhn+θhn-1,ϕh〉+𝒩(uhn+1,θhn+1,ϕh)+λμ〈∇θhn+1,∇ϕh〉=〈b(tn+1),ϕh〉,∀ϕh∈𝕋h.Remark 2 (discontinuity ofuhn+1).
We note thatuhn+1 is a discontinuous function across inter-element boundaries and that, in light of (7) and (8), uhn+1 is discrete divergence free in the sense that
(12)〈uhn+1,∇ϕh〉=〈u^hn+1+δ∇ψhn+1,∇ϕh〉=0,∀ϕh∈ℙh.We now summarize the results of this paper along with organization. We introduce appropriate Assumptions1–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:
(13)∥uhN+1∥02+∥u^hN+1∥02+∥2uhN+1-uhN∥02+∥2θhN+1-θhN∥02+∥θhN+1∥02+∑n=1N‍(∥δδuhn+1∥02+∥δδθhn+1∥02+∥∇δψhn+1∥02)+3∥∇ψhN+1∥02+2μτ∥qhN+1∥02+μτ∑n=1N‍(∥∇u^hn+1∥02+λ∥∇θhn+1∥)≤∥uh1∥02+∥2uh1-uh0∥02+∥2θh1-θh0∥02+∥θh1∥02+3∥∇ψh1∥02+2μτ∥qh1∥02+Cτ∑n=1N‍(∥f(tn+1)∥-12+∥b(tn+1)∥-12).We then will carry out the following optimal error estimates through several lemmas in Section4.Theorem 4 (error estimates).
Suppose the exact solution of (1) is smooth enough and τ≤Ch. If Assumptions 1, 3, 4, and 5 mentioned later hold, then the errors of Algorithm 1 will be bound by
(14)τ∑n=1N‍(∥u(tn+1)-uhn+1∥02+∥u(tn+1)-u^hn+1∥02)≤C(τ4+h4),τ∑n=1N‍(∥u(tn+1)-u^hn+1∥12+h2∥θ(tn+1)-θhn+1∥12)≤C(τ2+h2),τ∑n=1N‍∥θ(tn+1)-θhn+1∥02≤C(τ4+h4).
Moreover, if Assumption 2 also hold, then one has
(15)τ∑n=1N‍∥p(tn+1)-phn+1∥02≤C(τ2+h2).We note that the conditionτ≤Ch 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(16)-Δv+∇r=z,∇·v=0,inΩ,
which is the stationary Stokes system with vanishing boundary condition v=0, as well as Poisson’s equation
(17)-Δω=ξ,inΩ,
with boundary condition ω=0.We now state a basic assumption aboutΩ.Assumption 1 (regularity of(v,r) and ω).
The unique solutions{v,r} of (16) and ω of (17) satisfy
(18)∥v∥2+∥r∥1≤C∥z∥0,∥ω∥2≤C∥ξ∥0.
We remark that the validity of Assumption 1 is known if ∂Ω is of class C2 [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𝕍h and ℙh.Assumption 2 (discrete inf-sup).
There exists a constantβ>0 such that
(19)infsh∈ℙhsupwh∈𝕍h〈∇·wh,sh〉∥wh∥1∥sh∥0≥β.Assumption 3 (shape regularity and quasiuniformity [8–10]).
There exists a constantC>0 such that the ratio between the diameter hK of an element K∈𝔗 and the diameter of the largest ball contained in K is bounded uniformly by C, and hK is comparable with the meshsize h for all K∈𝔗.In order to launch Algorithm1, we need to set (uh1,ph1,θh1) via any first-order methods which hold the following conditions.Assumption 4 (the setting of the first step values).
Let(u(t1),p(t1),θ(t1)) be the exact solution of (1) at t=t1. The first step value (uh1,ph1,θh1) satisfies
(20)∥u(t1)-uh1∥0+∥θ(t1)-θh1∥0≤C(τ2+h2),∥u(t1)-uh1∥1+∥p(t1)-ph1∥1+∥θ(t1)-θh1∥1≤C(τ+h).Assumption 5 (approximability [8–10]).
For each(w,η,s)∈H2(Ω)×H2(Ω)×H1(Ω), there exist approximations (wh,ηh,sh)∈𝕍h×𝕋h×ℙh such that
(21)∥w-wh∥0+h∥w-wh∥1≤Ch2∥w∥2,∥η-ηh∥0+h∥η-ηh∥1≤Ch2∥η∥2,∥s-sh∥0≤Ch∥s∥1.The following elementary but crucial relations are derived in [14].Lemma 5 (inverse inequality).
IfIh denotes the Clément interpolant, then
(22)∥Ihw∥L3(Ω)≤Ch-d/6∥Ihw∥0,∥w-Ihw∥L3(Ω)≤Ch2-d/6∥w∥2.Lemma 6 (div-grad relation).
Ifw∈H01(Ω), then
(23)∥∇·w∥0≤∥∇w∥0.Let now(vh,rh)∈𝕍h×ℙh indicate the finite element solution of (16); namely,
(24)〈∇vh,∇wh〉+〈∇rh,wh〉=〈z,wh〉,∀wh∈𝕍h,〈∇·vh,sh〉=0,∀sh∈ℙh.Then we can find the well-known lemmas in [8–10].Lemma 7 (error estimates for mixed FEM).
Let(v,r)∈H01(Ω)×L02(Ω) be the solutions of (16), and, (vh,rh)=𝔖h(v,r)∈𝕍h×ℙh be the Stokes projection defined by (24), respectively. If Assumptions 1, 2, 3, and 5 hold, then
(25)∥v-vh∥0+h∥v-vh∥1+h∥r-rh∥0≤Ch2(∥v∥2+∥r∥1),(26)∥|v-vh|∥:=∥v-vh∥L∞(Ω)+∥∇(v-vh)∥L3(Ω)≤C∥z∥0.Proof.
Inequality (25) is standard [8–10]. To establish (26), we just deal with the L∞-norm, since the other can be treated similarly. If Ih denotes the Clément interpolant, then ∥v-Ihv∥L∞(Ω)≤C∥v∥2 and
(27)∥Ihv-vh∥L∞(Ω)≤Ch-d/2∥Ihv-vh∥L2(Ω)≤C∥v∥2
as a consequence of an inverse estimate. This completes the proof.Remark 8 (H1 stability of rh).
The bound∥∇rh∥0≤C(∥v∥2+∥r∥1) is a simple by-product of (16). To see this, we add and subtract Ihq, use the stability of Ih in H1, and observe that (25) implies ∥∇(rh-Ihr)∥0≤Ch-1∥rh-Ihr∥0≤C.Lemma 9 (error estimates for Poisson's equation [8, 10]).
Letω∈H01(Ω) be the solution of Poisson's equation (17), and, ωh∈𝕋h be its finite element approximation
(28)〈∇ωh,∇ϕh〉=〈ξ,ϕh〉,∀ϕh∈𝕋h.
If Assumptions 1, 3, and 5 hold, then there exists a positive constant C satisfying
(29)∥ω-ωh∥0+h∥ω-ωh∥1≤Ch2∥ω∥2≤Ch2∥ξ∥0,∥|ω-ωh|∥≤C∥ξ∥0.We finally state without proof several properties of the nonlinear form𝒩. In view of (5), we have the following properties of 𝒩 for all uh,vh,wh∈𝕍h:
(30)𝒩(uh,vh,wh)=-𝒩(uh,wh,vh),𝒩(uh,vh,vh)=0,(31)∇·u=0⇒𝒩(u,vh,wh)=〈(u·∇)vh,wh〉=-〈(u·∇)wh,vh〉.
Applying Sobolev imbedding lemma yields the following useful results.Lemma 10 (bounds on nonlinear convection [1, 12]).
Letu,v∈H2(Ω) with ∇·u=0, and let uh,vh,wh∈𝕍h. Then
(32)𝒩(u,vh,wh)≤C{∥u∥1∥vh∥1∥wh∥1∥u∥2∥vh∥1∥wh∥0∥u∥2∥vh∥0∥wh∥1,(33)𝒩(uh,v,wh)≤C∥uh∥0∥v∥2∥wh∥1.
In addition,
(34)𝒩(uh,vh,wh)≤C{∥uh∥0∥|vh|∥∥wh∥1∥uh∥L3(Ω)∥vh∥1∥wh∥1.Remark 11 (Treatment of convection term).
As we mention in Remark2, uhn+1 is discontinuous across inter-element boundaries and so uhn+1∉H1(Ω). 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{zn}n=0N, one has
(35)2〈3zn+1-4zn+zn-1,zn+1〉=δ∥zn+1∥02+δ∥2zn+1-zn∥02+∥δδzn+1∥02,(36)2〈zn+1-zn,zn+1〉=∥zn+1∥02-∥zn∥02+∥zn+1-zn∥02,(37)2〈zn+1-zn,zn〉=∥zn+1∥02-∥zn∥02-∥zn+1-zn∥02.


## 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:
(38)12τ〈3uhn+1-4uhn+uhn-1,wh〉+𝒩(uh*,u^hn+1,wh)+κμ2〈gθh*,wh〉-〈∇(32τψhn+1-μqhn),wh〉+μ〈∇u^hn+1,∇wh〉=〈f(tn+1),wh〉.
We now choose wh=4τu^hn+1 and apply (35) to obtain
(39)δ∥uhn+1∥02+δ∥2uhn+1-uhn∥02+∥δδuhn+1∥02+4τμ∥∇u^hn+1∥02=∑i=14‍Ai,
where
(40)A1:=6〈∇ψhn+1,u^hn+1〉,A2:=4τμ〈qhn,∇·u^hn+1〉,A3:=4τ〈f(tn+1),u^hn+1〉,A4:=-4τκμ2〈g(2θhn-θhn-1),u^hn+1〉.
We give thanks to (30) for eliminating the convection term. In light of u^hn+1=uhn+1-∇δψhn+1, (12) and (36) yield
(41)A1=-6〈∇ψhn+1,∇δψhn+1〉=-3(∥∇ψhn+1∥02-∥∇ψhn∥02+∥∇δψhn+1∥02).
Before we estimate A2, we evaluate an inequality via choosing ϕh=δqhn+1 in (9) to get
(42)∥δqhn+1∥02=-〈∇·u^hn+1,δqhn+1〉≤∥∇·u^hn+1∥0∥δqhn+1∥0.
Lemma 6 derives ∥δqhn+1∥02≤∥∇·u^hn+1∥02≤∥∇u^hn+1∥02, and so (9) and (37) lead us to
(43)A2=-4μτ〈qhn,δqhn+1〉=-2μτ(∥qhn+1∥02-∥qhn∥02-∥δqhn+1∥02)≤-2μτ(∥qhn+1∥02-∥qhn∥02)+2μτ∥∇u^hn+1∥02.
Clearly, we have
(44)A3≤Cτμ∥f(tn+1)∥-12+τμ∥∇u^hn+1∥02.
We now attack A4 with ∥u^hn+1∥02=∥uhn+1∥02+∥∇δψhn+1∥02 which comes from (8) and (12). Then we arrive at
(45)A4≤Cκ2μ4τ∥2θhn-θhn-1∥02+Cτ∥uhn+1∥02+∥∇δψhn+1∥02.
In conjunction with 𝒩(uh*,θhn+1,θhn+1)=0 which comes from (30), if we choose ϕh=4τθhn+1 in (11) and use (35), then we obtain
(46)δ∥θhn+1∥02+δ∥2θhn+1-θhn∥02+∥δδθhn+1∥02+4τλμ∥∇θhn+1∥=4τ〈b(tn+1),θhn〉≤Cτλμ∥b(tn+1)∥-12+λμτ∥∇θhn∥02.
Inserting A1–A4 back into (39), adding with (46), and then summing over n from 1 to N lead to (13) by help of discrete Gronwall inequality and the equality ∥u^hn+1∥02=∥uhn+1∥02+∥∇δψhn+1∥02. So we finish the proof of Theorem 3.


## 4. Error Estimates

We prove Theorem4 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 (Uhn+1,Phn+1):=𝔖h(u(tn+1),p(tn+1))∈𝕍h×ℙh to be the Stokes projection of the true solution at time tn+1. It means that (Uhn+1,Phn+1)∈𝕍h×ℙh is the solution of
(47)〈∇(u(tn+1)-Uhn+1),∇wh〉+〈∇(p(tn+1)-Phn+1),wh〉=0,∀wh∈𝕍h,〈∇·Uhn+1,ϕh〉=0,ϕh∈ℙh.
We also define Θhn+1∈𝕋h as the solution of
(48)〈∇Θhn+1,∇ϕh〉=〈∇θ(tn+1),∇ϕh〉,∀ϕh∈𝕋h.
And we denote notations Gn+1:=u(tn+1)-Uhn+1,
(49)gn+1:=p(tn+1)-Phn+1,ηn+1:=θ(tn+1)-Θhn+1.
From Lemma 7, we can deduce
(50)∥Gn+1∥02+h2∥Gn+1∥12+h2∥gn+1∥02≤Ch4(∥u(tn+1)∥22+∥p(tn+1)∥12),(51)∥δGn+1∥02+h2∥δGn+1∥12+h2∥δgn+1∥02≤Cτh4∫tntn+1‍(∥ut(t)∥22+∥pt(t)∥12)dt.
In conjunction with the definition of |∥·∥| in (26), we can derive
(52)∥|Gn+1|∥≤C.
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
(53)Ehn+1:=Uhn+1-uhn+1,E^hn+1:=Uhn+1-u^hn+1,ehn+1:=Phn+1-phn+1,ϑhn+1:=Θhn+1-θhn+1.
Then, in conjunction with (50), we can readily get the same accuracy for the errors
(54)En+1:=u(tn+1)-uhn+1,E^n+1:=u(tn+1)-u^hn+1,en+1:=p(tn+1)-phn+1,ϑn+1:=θ(tn+1)-θhn+1.
Additionally, we denote
(55)εhn+1:=Phn+1+3ψhn+12τ.
We readily obtain the following properties:
(56)Gn+1=En+1-Ehn+1=E^n+1-E^hn+1,E^hn+1=Ehn+1+∇δψhn+1,E^n+1=E^hn+1=Gn+1=0,on∂Ω.
Moreover, from (12),
(57)〈Ehn+1,∇ϕh〉=〈En+1,∇ϕh〉=〈Gn+1,∇ϕh〉=0,∀ϕh∈ℙh.
Whence we deduce crucial orthogonality properties:
(58)∥E^hn+1∥02=∥Ehn+1∥02+∥∇δψhn+1∥02.
We also point out that, owing to Lemma 6, qhn+1∈ℙh defined in (9) satisfies
(59)∥qhn+1-qhn∥0≤∥∇E^n+1∥0.
In conjunction with ηn+1=θ(tn+1)-Θhn+1=ϑn+1-ϑhn+1, Assumption 5 leads to
(60)∥ηn+1∥02+h2∥ηn+1∥1≤Ch4∥θ(tn+1)∥22,∥δηn+1∥02+h2∥δηn+1∥12≤Cτh4∫tntn+1‍∥θt∥2dt.We now estimate the first-order accuracy for velocity and temperature in Lemma13 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
(61)∥E^hN+1∥02+∥EhN+1∥02+∥2EhN+1-EhN∥02+2μτ∥qhN+1∥02+4τ23∥∇εhN+1∥02+∥ϑhN+1∥02+∥2ϑhN+1-ϑhN∥02+μτ∑n=1N‍(∥∇E^hn+1∥02+λ∥∇ϑhn+1∥02)+∑n=1N‍(∥δδEhn+1∥02+∥∇δψhn+1∥02+∥δδϑhn+1∥02)≤C(τ2+h2).Proof.
We resort to the Taylor theorem to write (1) as follows:
(62)3u(tn+1)-4u(tn)+u(tn-1)2τ+(u(tn+1)·∇)u(tn+1)+∇p(tn+1)-μΔu(tn+1)+κμ2gθ(tn+1)=Rn+1+f(tn+1),(63)3θ(tn+1)-4θ(tn)+θ(tn-1)2τ+u(tn+1)·∇θ(tn+1)-λμΔθ(tn+1)=Qn+1+b(tn+1),
where
(64)Rn+1:=(1τ)∫tntn+1‍uttt(s)(s-tn)2ds-(14τ)∫tn-1tn+1‍uttt(s)(tn-1-s)2ds,Qn+1:=(1τ)∫tntn+1‍θttt(s)(s-tn)2ds-(14τ)∫tn-1tn+1‍θttt(s)(tn-1-s)2ds
are the truncation errors. In conjunction with the definition of the Stokes projection {Uhn+1,Phn+1}, we readily get a weak formulation for (62), ∀wh∈𝕍h:
(65)12τ〈3u(tn+1)-4u(tn)+u(tn-1),wh〉+𝒩(u(tn+1),u(tn+1),wh)+〈∇Phn+1,wh〉+μ〈∇Uhn+1,∇wh〉+κμ2〈gθ(tn+1),wh〉=〈Rn+1,wh〉+〈f(tn+1),wh〉.
We replace the pressure phn term in (6) with (10) and then subtract it from (65) to obtain
(66)12τ〈3E^n+1-4En+En-1,wh〉+μτ〈∇E^hn+1,∇wh〉=-〈∇(δPhn+1+εhn),wh〉+μ〈∇qhn,wh〉+𝒩(2uhn-uhn-1,u^hn+1,wh)-𝒩(u(tn+1),u(tn+1),wh)-κμ2〈g(θ(tn+1)-2θhn+θhn-1),wh〉+〈Rn+1,wh〉.
Choosing wh=4τE^hn+1=4τ(E^n+1-Gn+1)=4τ(Ehn+1+∇δψhn+1) in (66) and using (58) and (35) yield
(67)δ∥Ehn+1∥02+δ∥2Ehn+1-Ehn∥02+∥δδEhn+1∥02+6∥∇δψhn+1∥02+4μτ∥∇E^hn+1∥02=∑n=17‍Ai,
where
(68)A1:=-4τ𝒩(u(tn+1),u(tn+1),E^hn+1)+4τ𝒩(2uhn-uhn-1,u^hn+1,E^hn+1),A2:=-2〈3Gn+1-4Gn+Gn-1,E^hn+1〉,A3:=-4τ〈∇δPhn+1,E^hn+1〉,A4:=-4τ〈∇εhn,E^hn+1〉,A5:=4μτ〈∇qhn,E^hn+1〉,A6:=4τ〈Rn+1,E^hn+1〉,A7:=-4τκμ2〈g(θ(tn+1)-2θhn+θhn-1),E^hn+1〉.
We now estimate all the terms from A1 to A7, respectively.  To tackle A1, we first add and subtract 2u(tn)-u(tn-1) to obtain
(69)A1=-4τ𝒩(δδu(tn+1),u(tn+1),E^hn+1)-4τ𝒩(2uhn-uhn-1,E^n+1,E^hn+1)-4τ𝒩(2En-En-1,u(tn+1),E^hn+1).
Because of 𝒩(2uhn-uhn-1,E^hn+1,E^hn+1)=0, which comes from (30), the second term of A1 can be replaced by
(70)4τ𝒩(2En-En-1-2u(tn)+u(tn-1),Gn+1,E^hn+1).
If we apply Lemma 10, then we can readily obtain
(71)A1≤Cτ(∥δδu(tn+1)∥0∥u(tn+1)∥2+∥2En-En-1∥0∥u(tn+1)∥2)∥E^hn+1∥1+Cτ(∥2En-En-1∥0∥|Gn+1|∥+∥2u(tn)-u(tn-1)∥2∥Gn+1∥0∥|Gn+1|∥)∥E^hn+1∥1.
Since we have ∥u(tn+1)∥2+∥|Gn+1|∥≤M according to (52), we arrive at
(72)A1≤Cτ(∥2Ehn-Ehn-1∥02+∥2Gn-Gn-1∥02+∥Gn+1∥02∥2Ehn-Ehn-1∥02)+μτ∥∇E^hn+1∥02+Cτ4μ×∫tn-1tn+1‍∥uttt(t)∥02dt.
In light of (51) and (58), A2 becomes
(73)A2≤Cτ∥Ehn+1∥02+14∥∇δψhn+1∥02+Ch4∫tn-1tn+1‍(∥ut(t)∥22+∥pt(t)∥12)dt.
In order to estimate A3 and A4, we note that the cancellation law (12) gives 〈∇εhn,Ehn+1〉=0. Then E^hn+1=Ehn+1+∇δψhn+1 yields
(74)A3=-4τ〈∇δPhn+1,∇δψhn+1〉≤14∥∇δψhn+1∥02+Cτ3×∫tntn+1‍(∥ut(t)∥22+∥pt(t)∥12)dt.
In conjunction with the definition εhn+1=Phn+1+3ψhn+1/2τ, A4 can be evaluated by
(75)A4=-4τ〈∇εhn,∇δψhn+1〉=-8τ23〈∇εhn,∇(δεhn+1-δPhn+1)〉≤-4τ23(∥∇εhn+1∥02-∥∇εhn∥02-∥∇δεhn+1∥02)+Cτ3∥∇εhn∥02+Cτ∥∇δPhn+1∥02.
If we now apply inequality (a+b)2≤4a2+(4/3)b2, then we can get
(76)4τ23∥∇δεhn+1∥02=4τ23∥∇δPhn+1+32τ∇δψhn+1∥02≤Cτ2∥∇δPhn+1∥02+4∥∇δψhn+1∥02.
So we arrive at
(77)A4≤-4τ23(∥∇εhn+1∥02-∥∇εhn∥02)+Cτ3∥∇εhn∥02+4∥∇δψhn+1∥02+Cτ2∫tntn+1‍(∥ut(t)∥22+∥pt(t)∥12)dt.
In light of ∇·u(tn+1)=0 and (37), (9) and (59) yield
(78)A5=4μτ〈qhn,∇·u^hn+1〉=-4μτ〈qhn,δqhn+1〉=-2μτ(∥qhn+1∥02-∥qhn∥02-∥δqhn+1∥02)≤-2μτ(∥qhn+1∥02-∥qhn∥02)+2μτ∥∇E^hn+1∥02.
Also we readily get
(79)A6≤Cτ∥Ehn+1∥02+14∥∇δψhn+1∥02+Cτ4∫tn-1tn+1‍∥uttt(t)∥02dt.
The Hölder inequality and (60) yield
(80)A7=-4κμ2τ〈g(δδθ(tn+1)-2ϑn+ϑn-1),E^hn+1〉≤Cκ2μ4τ(∥2ϑhn-ϑhn-1∥02+∥2ηn-ηn-1∥02)+τ∥Ehn+1∥02+14∥∇δψhn+1∥02+Cτ4∫tn-1tn+1‍∥θtt(t)∥02dt.
Inserting the above estimates into (67) gives
(81)δ∥Ehn+1∥02+δ∥2Ehn+1-Ehn∥02+∥δδEhn+1∥02+∥∇δψhn+1∥02+μτ∥∇E^hn+1∥02+4τ23∥∇εhn+1∥02+2μτ∥qhn+1∥02≤4τ23∥∇εhn∥02+Cτ3∥∇εhn∥02+2μτ∥qhn∥02+Cτ(∥Ehn+1∥02+∥2Ehn-Ehn-1∥02+∥2Gn-Gn-1∥02+∥Gn+1∥02∥Ehn+1∥02)+Cκ2μ4τ∥2ηn-ηn-1∥02+Cτ4∫tn-1tn+1‍(∥uttt(t)∥02+∥θtt(t)∥02)dt+Cκ2μ4τ∥2ϑhn-ϑhn-1∥02+C(τ2+h4)∫tn-1tn+1‍(∥ut(t)∥22+∥pt(t)∥12)dt.
On the other hand, the definition (48) of Θhn+1∈𝕋h leads to a weak formula of (63) as, for all ϕh∈𝕋h,
(82)12τ〈3θ(tn+1)-4θ(tn)+θ(tn-1),ϕh〉+𝒩(u(tn+1),θ(tn+1),ϕh)+λμ〈∇Θhn+1,∇ϕh〉=〈Qn+1+b(tn+1),ϕh〉.
We now subtract (11) from (82) to derive
(83)12τ〈3ϑn+1-4ϑn+ϑn-1,ϕh〉+𝒩(u(tn+1),θ(tn+1),ϕh)-𝒩(uhn+1,θhn+1,ϕh)+λμ〈∇ϑhn+1,∇ϕh〉=〈Qn+1,ϕh〉.
Choosing ϕh=4τϑhn+1=4τ(ϑn+1-ηn+1) yields
(84)δ∥ϑhn+1∥02+δ∥2ϑhn+1-ϑhn∥02+∥δδϑhn+1∥02+4λμτ∥∇ϑhn+1∥02=A8+A9,
where
(85)A8:=-4τ𝒩(u(tn+1),θ(tn+1),ϑhn+1)+4τ𝒩(uhn+1,θhn+1,ϑhn+1),A9:=-2〈3ηn+1-4ηn+ηn-1,ϑhn+1〉+4τ〈Qn+1,ϑhn+1〉.
To estimate A8, we note 𝒩(En+1-u(tn+1),ϑhn+1,ϑhn+1)=0 which comes from (30), then ∥En+1∥L3(Ω)≤Ch-1/2∥En+1∥0 and ∥ηn+1∥1≤Ch∥θ(tn+1)∥2 yield
(86)A8=-4τ𝒩(En+1,θ(tn+1),ϑhn+1)+4τ𝒩(En+1-u(tn+1),ηn+1,ϑhn+1)≤Cτ∥En+1∥0∥θ(tn+1)∥2∥ϑhn+1∥1+Cτ(∥En+1∥L3(Ω)∥ηn+1∥1+∥u(tn+1)∥2∥ηn+1∥0En+1)∥ϑhn+1∥1≤Cτλμ(∥En+1∥02+∥ηn+1∥02)+λμτ∥∇ϑhn+1∥02.
We use ∥δηn+1∥02≤Cτh4∫tntn+1‍∥θt(t)∥22dt which is (60) to attack A9:
(87)A9≤C∥3ηn+1-4ηn+ηn-1∥0∥ϑhn+1∥0+Cτ∥Qn+1∥0∥ϑhn+1∥0≤Cτ∥ϑhn+1∥02+Ch4∫tn-1tn+1‍∥θt(t)∥22dt+Cτ4∫tn-1tn+1‍∥θttt(t)∥02dt.
Inserting the above estimates into (84) gives
(88)δ∥ϑhn+1∥02+δ∥2ϑhn+1-ϑhn∥02+∥δδϑhn+1∥02+3λμτ∥∇ϑhn+1∥02≤Cτ∥ϑhn+1∥02+Cτλμ(∥En+1∥02+∥ηn+1∥02)+C(h4+τ4)∫tn-1tn+1‍(∥θttt(t)∥02+∥θt(t)∥22)dt.
Adding 2 inequalities (81) and (88) and summing over n from 1 to N imply
(89)∥EhN+1∥02+∥2EhN+1-EhN∥02+∥ϑhN+1∥02+∥2ϑhN+1-ϑhN∥02+2μτ∥qhn+1∥02+4τ23∥∇εhn+1∥02+∑n=1N‍(∥δδEhN+1∥02+∥δδϑhn+1∥02+∥∇δψhn+1∥02)+μτ∑n=1N‍(∥∇E^hn+1∥02+3λ∥∇ϑhn+1∥02)≤∥Eh1∥02+∥2Eh1-Eh0∥02+∥ϑh1∥02+∥2ϑh1-ϑh0∥02+4τ23∥∇εh1∥02+2μτ∥qh1∥02+τ∑n=1N‍(∥En+1∥02+∥ηn+1∥02)+Cτ∑n=1N‍(∥ϑhn+1∥02+∥2ϑhn-ϑhn-1∥02∥ϑhn+1∥02+∥2ηn-ηn-1∥02)+Cτ3∥∇εhn∥02+C(h4+τ4)×∫tn-1tn+1‍(∥uttt(t)∥02+∥θttt(t)∥02+∥θt(t)∥22)dt+C(τ2+h4)∫tn-1tn+1‍(∥ut(t)∥22+∥pt(t)∥12)dt.
Because of q1=0, we have εh1=Ph1+(3/2τ)ψh1=eh1, and thus, Assumption 4 yields ∥∇εh1∥02≤C and the first 4 terms in the right-hand side can be bound by Assumption 4 and properties E0=0 and q1=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 ofA3 and A4 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 n=1, and so we have to recompute again (77). We start to rewrite A4 as
(90)A4=-8τ23〈∇εh1,∇(δεh2-δPh2)〉≤-4τ23(∥∇εh2∥02-∥∇εh1∥02-∥∇δεh2∥02)+Cτ2∥∇εh1∥02+Cτ2∥∇δPhn+1∥02≤-4τ23∥∇εh2∥02+Cτ2∥∇εh1∥02+4∥∇δψh2∥02+Cτ3×∫t0t2(∥ptt(t)∥12+∥utt(t)∥22)‍dt.
In light of Assumption 4, we arrive at
(91)∥E^h2∥02+∥Eh2∥02+∥2Eh2-Eh1∥02+12∥δδEh2∥02+∥∇δψh2∥02+∥2ϑh2-ϑh2∥02+∥ϑh2∥02+μτ∥∇E^h2∥02+2μτ∥qh2∥02+4τ23∥∇εh2∥02+μλτ∥∇ϑh2∥02≤C(τ4+h4).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 τ≤Ch. If Assumptions 3–5 hold, then the time derivative velocity and temperature error functions satisfy
(92)∥δEhN+1∥02+∥δE^hN+1∥02+∥2δEhN+1-δEhN∥02+2μτ∥δqhN+1∥02+∥δϑhN+1∥02+∥2δϑhN+1-δϑhN∥02+∑n=2N(∥δδδEhn+1∥02+∥∇δδψhn+1∥02‍+∥δδδϑhn+1∥02)+4τ23∥∇δεhN+1∥02+μτ∑n=2N‍(∥∇δE^hn+1∥02+λ∥∇δϑhn+1∥02)≤C(τ2+h2)2.Remark 16 (the conditionτ≤Ch).
The assumptionτ≤Ch 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 h.Proof.
We subtract two consecutive formulas of (66) and impose wh=4τδE^hn+1 to obtain
(93)∥δEhn+1∥02+∥2δEhn+1-δEhn∥02+∥δδδEhn+1∥02+6∥∇δδψhn+1∥02-∥δEhn∥02-∥2δEhn-δEhn-1∥02+4μτ∥∇δE^hn+1∥02=∑n=17‍Ai,
where
(94)A1:=4τ𝒩(u(tn),u(tn),δE^hn+1)-4τ𝒩(2uhn-1-uhn-2,u^hn,δE^hn+1)-4τ𝒩(u(tn+1),u(tn+1),δE^hn+1)+4τ𝒩(2uhn-uhn-1,u^hn+1,δE^hn+1),A2:=-2〈3δGn+1-4δGn+δGn-1,δE^hn+1〉,A3:=-4τ〈∇δδPhn+1,δE^hn+1〉,A4:=-4τ〈∇δεhn,δE^hn+1〉,A5:=4μτ〈∇δqhn,δE^hn+1〉,A6:=4τ〈δRn+1,δE^hn+1〉,A7:=4κμ2τ×〈gδ(θ(tn+1)-2θn+θn-1),δE^hn+1〉.
We now estimate each term from A1 to A7  separately. The convection term A1 can be rewritten as follows:
(95)A1=4τ𝒩(δδu(tn),u(tn),δE^n+1)-4τ𝒩(δδu(tn+1),u(tn+1),δE^n+1)+4τ𝒩(2En-1-En-2,u(tn),δE^n+1)-4τ𝒩(2En-En-1,u(tn+1),δE^n+1)+4τ𝒩(2uhn-1-uhn-2,E^n,δE^n+1)-4τ𝒩(2uhn-uhn-1,E^n+1,δE^n+1).
In estimating convection terms, we will use Lemma 10 frequently without notice. We recall ∥u(t)∥2≤C to obtain
(96)A1,1+A1,2≤Cτ(∥δδu(tn+1)∥0∥u(tn+1)∥2+∥δδu(tn)∥0∥u(tn)∥2tn+1)∥δE^hn+1∥1≤μτ6∥∇δE^hn+1∥02+Cτ4μ∫tn-2tn+1‍∥utt(t)∥02dt.
The result in Lemma 13, ∥2En-En-1∥0≤C(τ+h), is essential to treat the next 2 convection terms. Invoking (33), we have
(97)A1,3+A1,4≤Cτ∥2En-En-1∥0×∥δu(tn+1)∥2∥δE^hn+1∥1+Cτ∥2δEn-δEn-1∥0×∥u(tn)∥2∥δE^hn+1∥1≤μτ6∥∇δE^hn+1∥02+Cτμ∥2δEn-δEn-1∥02+Cτ2(τ2+h2)μ×∫tn-1n‍∥ut(t)∥22dt.
We note 𝒩(2uhn-uhn-1,δE^hn+1,δE^hn+1)=0 which comes from (30). Then we obtain
(98)A1,5+A1,6=-4τ𝒩(E^hn+12δuhn-δuhn-1,Gn+1+E^hn+1,δE^hn+1)-4τ𝒩(2uhn-1-uhn-2,δGn+1,δE^hn+1)=4τ𝒩(2δEn-δEn-1-2δu(tn)+δu(tn-1),Gn+1+E^hn+1,δE^hn+1)+4τ𝒩(2En-1-En-2-2u(tn-1)+u(tn-2),δGn+1,δE^hn+1):=B1+B2.
To attack B1, we first note Lemma 5 which is, for any wh∈𝕍h, ∥wh∥L3(Ω)≤Ch-d/6∥wh∥0. If we apply
(99)∥E^n+1+Gn+1∥0+τ+h∥E^n+1+Gn+1∥1≤C(τ+h)
which is the result of Lemma 13, then we can conclude, in light of (34),
(100)B1≤Cτ∥2δEn-δEn-1∥L3(Ω)×∥E^hn+1+Gn+1∥1∥∇δE^hn+1∥0+Cτ∥2δu(tn)-δu(tn-1)∥2×∥E^hn+1+Gn+1∥0∥∇δE^hn+1∥0≤μτ6∥∇δE^hn+1∥02+Cτμ(1+τh)×∥2δEn-δEn-1∥02+Cτ2μ(τ+h)2∫tn-2tn‍∥ut(t)∥22dt.
We now estimate B2 using ∥2En-1-En-2∥L3(Ω)≤Ch-d/6∥2En-1-En-2∥0≤M:
(101)B2≤Cτ∥2En-1-En-2∥L3(Ω)×∥δGn+1∥1∥∇δE^hn+1∥0+Cτ∥2u(tn-1)-u(tn-2)∥1×∥δGn+1∥1∥∇δE^hn+1∥0≤Cτ∥δGn+1∥12+μτ6∥∇δE^hn+1∥02≤μτ6∥∇δE^hn+1∥02+Cτ2h2×∫tn-1tn+1‍(∥ut(t)∥22+∥pt(t)∥12)dt.
In light of Hölder inequality, (51) yields
(102)A2=-2〈3δGn+1-4δGn+δGn-1,δE^hn+1〉≤μτ6∥∇δEhn+1∥02+Ch4∫tn-2tn+1‍(∥ut(t)∥22+∥pt(t)∥12)dt.
Integral by parts leads to
(103)A3=4τ〈δδPhn+1,∇·δE^hn+1〉≤μτ6∥∇δE^hn+1∥02+Cτ4μ×∫tn-1tn+1‍(∥utt(t)∥22+∥ptt(t)∥12)dt.
In order to tackle A4, marking use of δE^hn+1=δEhn+1+∇δδψhn+1, we readily get
(104)A4=-4τ〈∇δεhn,δEhn+1+∇δδψhn+1〉=-8τ23〈∇δεhn,∇(δδεhn+1-δδPhn+1)〉≤-4τ23(∥∇δεhn+1∥02-∥∇δεhn∥02-∥∇δδεhn+1∥02)+Cτ3∥∇δεhn∥02+Cτ∥∇δδPhn+1∥02.
If we now apply inequality (a+b)2≤4a2+(4/3)b2, then we can have
(105)4τ23∥∇δδεhn+1∥02=4τ23∥∇δδPhn+1+32τ∇δδψhn+1∥02≤Cτ2∥∇δδPhn+1∥02+4∥∇δδψhn+1∥02.
So we arrive at
(106)A4≤-4τ23(∥∇δεhn+1∥02-∥∇δεhn∥02)+4∥∇δδψhn+1∥02+Cτ3∥∇δεhn∥02+Cτ4∫tn-1tn+1‍(∥utt(t)∥22+∥ptt(t)∥12)dt.
Invoking (9), (37) and (60) lead to
(107)A5=-4μτ〈δqhn,δδqhn+1〉=-2μτ(∥δqhn+1∥02-∥δqhn∥02-∥δδqhn+1∥02)≤-2μτ(∥δqhn+1∥02-∥δqhn∥02)+2μτ∥∇δE^hn+1∥02.
Finally, the A6 term becomes
(108)A6≤Cτ∥δRn+1∥0∥δE^hn+1∥0≤Cτ∥δE^hn+1∥02+Cτ4∫tn-1tn+1‍∥uttt(t)∥02dt.
For the new term A7, we argue as in Lemma 13 to arrive at
(109)A7≤Cκ2μ4τ×(∥2δϑhn-δϑhn-1∥02+∥2δηn-δηn-1∥02ϑhn-1)+τ∥δEhn+1∥02+12∥∇δδψhn+1∥02+Cτ4∫tn-2tn+1‍∥θtt(t)∥02dt.
We now insert the above estimates into (93) to obtain
(110)δ∥δEhn+1∥02+δ∥2δEhn+1-δEhn∥02+∥δδδEhn+1∥02+μτ∥∇δE^hn+1∥02+∥∇δδψhn+1∥02+2μτ∥δqhn+1∥02+4τ23∥∇δεhn+1∥02≤Cτμ(∥2δEhn-δEhn-1∥02Ehn-1+∥2δGn-δGn-1∥02)+4τ23∥∇δεhn∥02+Cκ2μ4τ×(∥2δϑhn-δϑhn-1∥02ϑhn-1+∥2δηn-δηn-1∥02)+τ∥δEhn+1∥02+2μτ∥δqhn∥02+Cτ4×∫tn-1tn+1(∥uttt(t)∥02+∥utt(t)∥22‍+∥ptt(t)∥12+∥θtt(t)∥02)dt+Cτ3∥∇δεhn∥02+C(τ2+h2)2×∫tn-2tn+1‍(∥utt(t)∥02+∥ut(t)∥22+∥pt(t)∥12)dt.
To evaluate errors ofδϑhn+1, we subtract two consecutive formulas (83) and choose by ϕh=4τδϑhn+1=4τ(δϑn+1-δηn+1) as follows:
(111)δ∥δϑhn+1∥02+δ∥2δϑhn+1∥02+∥δδδϑhn+1∥02+4λμτ∥∇δϑhn+1∥02=A8+A9,
where
(112)A8:=4τ𝒩(u(tn),θ(tn),δϑhn+1)-4τ𝒩(uhn,θhn,δϑhn+1)-4τ𝒩(u(tn+1),θ(tn+1),δϑhn+1)+4τ𝒩(uhn+1,θhn+1,δϑhn+1),A9:=-2〈3δηn+1-4δηn+δηn-1,δϑhn+1〉+4τ〈δQn+1,δϑhn+1〉.
We treat the A8 term first by rewriting it as follows:
(113)A8=-4τ𝒩(En+1,θ(tn+1),δϑhn+1)-4τ𝒩(uhn+1,ϑn+1,δϑhn+1)+4τ𝒩(En,θ(tn),δϑhn+1)+4τ𝒩(uhn,ϑn,δϑhn+1)∶=∑i=14‍A8,i.
In estimating convection terms, we will use Lemma 10 frequently without notice. We recall ∥θ(t)∥2≤C and ∥En∥0≤C(τ+h) which is the result in Lemma 13 to obtain
(114)A8,1+A8,3=-4τ𝒩(δEn+1,θ(tn+1),δϑhn+1)-4τ𝒩(En,δθ(tn+1),δϑhn+1)≤Cτ(∥δEn+1∥0∥θ(tn+1)∥2+∥En∥0∥δθ(tn+1)∥2)∥δϑhn+1∥1≤λμτ∥∇δϑhn+1∥02+Cτλμ∥δEn+1∥02+Cτ2(τ+h)2λμ×∫tntn+1‍∥θt(t)∥22dt.
In conjunction with 𝒩(uhn+1,δϑhn+1,δϑhn+1)=0 which comes from (30), we rewrite A8,2+A8,4 as
(115)A8,2+A8,4=4τ𝒩(Ehn+1-Uhn+1,δηn+1,δϑhn+1)+4τ𝒩(δEhn+1-δUhn+1,ϑn,δϑhn+1):=B3+B4.
Because we can readily get ∥Ehn+1-Uhn+1∥L3(Ω)≤M via Lemmas 5 and 13, we conclude, in conjunction with (60),
(116)B3≤Cτ∥Ehn+1-Uhn+1∥L3(Ω)×∥δηn+1∥1∥δϑhn+1∥1≤λμτ∥∇δϑhn+1∥02+Cτ2h2λμ×∫tntn+1‍∥θt(t)∥22dt.
Before we attack B4, we note that the assumption τ≤Ch is required to apply ∥δEhn∥L3(Ω)≤Cτ-1/2∥δEhn∥0. We now imply ∥ϑn∥02+τ∥ϑn∥12≤C(τ2+h2) and (51) to get
(117)B4=4τ𝒩(δGn+1-δu(tn+1)+δEhn+1,ϑn,δϑhn+1)≤Cτ((∥δGn+1∥1+∥δEhn+1∥L3(Ω))∥δEhn+1∥L3(Ω)×∥ϑn∥1+∥δu(tn)∥2∥ϑn∥0)∥δϑhn+1∥1≤Cτλμ∥δEhn+1∥02+λμτ∥∇δϑhn+1∥02+Cτ2(τ2+h2)λμ×∫tn-1tn+1(∥ut(t)∥22+∥pt(t)∥12)dt.‍
In light of Hölder inequality, A9 becomes
(118)A9≤Cτ∥δϑhn+1∥02+C(τ4+h4)×∫tn-2tn+1‍(∥θt(t)∥22+∥θttt∥02)dt.
Inserting the above estimates into (111) yields
(119)δ∥δϑhn+1∥02+δ∥2δϑhn+1-δϑhn∥02+∥δδδϑhn+1∥02+λμτ∥∇δϑhn+1∥02≤Cτ∥δϑhn+1∥02+Cτλμ(∥δEhn+1∥0+∥δEhn∥02+∥δGn+1∥0)+C(τ2+h2)2×∫tn-2tn+1‍(∥ut(t)∥22+∥pt(t)∥12+∥θt(t)∥22+∥θttt∥02)dt.
Adding 2 equations (110) and (119) and then summing up n from 2 to N 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(vn+1,rn+1) and (vhn+1,rhn+1) are the solutions of (16) and (24) with z=Ehn+1, respectively. And let ωn+1 and ωhn+1 be solutions of (17) and (28) with ξ=ϑhn+1. Let the exact solution of (1) be smooth enough and τ≤Ch. If Assumptions 1 and 3–5 hold, then one has
(120)∥∇vhN+1∥02+∥∇(2vhN+1-vhN)∥02+∥∇ωhN+1∥02+∥∇(2ωhN+1-ωhN)∥02+∑n=1N‍(∥∇δδvhn+1∥02+∥∇δδωhn+1∥02)+2μτ∑n=1N‍(∥Ehn+1∥02+λ∥ϑhn+1∥02)≤C(τ4+h4).Proof.
We choosewh=4τvhn+1∈𝕍h in (66) to obtain
(121)δ∥∇vhn+1∥02+δ∥∇(2vhn+1-vhn)∥02+∥∇δδvhn+1∥02+4μτ∥Ehn+1∥02=∑i=15‍Ai,
where
(122)A1:=4τ𝒩(2uhn-uhn-1,u^hn+1,vhn+1)-4τ𝒩(u(tn+1),u(tn+1),vhn+1),A2:=〈3Gn+1-4Gn+Gn-1,vhn+1〉,A3:=4μτ〈E^hn+1,∇rhn+1〉,A4:=4τ〈Rn+1,vhn+1〉,A5:=4κμ2τ〈g(θ(tn+1)-2θhn+θhn-1),vhn+1〉.
We now estimate all the terms from A1 to A5, respectively. The convection term A1 can be rewritten as follows:
(123)A1=-4τ𝒩(δδu(tn+1),u(tn+1),vhn+1)-4τ𝒩(2u(tn)-u(tn-1),E^n+1,vhn+1)+4τ𝒩(2En-En-1,E^n+1,vhn+1)-4τ𝒩(2En-En-1,u(tn+1),vhn+1),
and we denote by A1,i, for i=1,2,…,4, the four terms in the right-hand side. To estimate convection terms, we will use frequently Lemma 10 without notice. Using ∥u(tn+1)∥2≤M, we can readily get
(124)A1,1≤Cτ∥δδu(tn+1)∥0×∥u(tn+1)∥2∥vhn+1∥1≤Cτ∥∇vhn+1∥02+Cμτ4×∫tn-1tn+1‍∥utt(t)∥02dt,(125)A1,4≤Cτ(∥En∥0+∥δEn∥0)×∥u(tn+1)∥2∥vhn+1∥1≤μτ2(∥Ehn∥02+∥Gn∥02+∥δEhn∥02+∥δGn∥02)+Cτμ∥∇vhn+1∥02.
Because ∇·(2u(tn)-u(tn-1))=0 and 2u(tn)-u(tn-1)=0 on boundary, we can use (32) and so we get
(126)A1,2≤Cτ∥2u(tn)-u(tn-1)∥2×∥E^n+1∥0∥vhn+1∥1≤μτ2(∥Ehn+1∥02+∥Gn+1∥02+∥∇δψhn+1∥02)+Cτμ∥∇vhn+1∥02.
In light of εhn+1=Phn+1+3ψhn+1/2τ, we can obtain
(127)∥∇δψhn+1∥02=4τ29∥∇(δεhn+1-δPhn+1)∥02≤Cτ2∥∇δεhn+1∥02+Cτ3∫tntn+1‍(∥ut(t)∥22+∥pt(t)∥12)dt,
and so we can conclude
(128)A1,2≤μτ2(∥Ehn+1∥02+∥Gn+1∥02)+Cμτ3∥∇δεhn+1∥02+Cμτ4×∫tntn+1‍(∥ut(t)∥22+∥pt(t)∥12)dt+Cτμ∥∇vhn+1∥02.
If we apply ∥2En-En-1∥L3(Ω)≤Ch-d/6∥2En-En-1∥0≤Ch-d/6(τ+h) which can be derived from Lemmas 5 and 13, then we can get, by the help of Lemma 10 and Assumption 1,
(129)A1,3≤Cτ∥2En-En-1∥0×∥E^n+1∥1∥vn+1∥2+Cτ∥2En-En-1∥L3(Ω)×∥E^n+1∥1∥vn+1-vhn+1∥1≤Cτ(τ+h)∥E^n+1∥1∥vn+1∥2+Cτh∥E^n+1∥1∥vn+1∥2≤Cμτ(τ+h)2×(∥∇E^hn+1∥02+∥∇Gn+1∥02)+μτ2∥Ehn+1∥02.
In conjunction with (51), we can have
(130)A2≤Cτ(∥δGn+1∥02+∥δGn∥02)+Cτ∥vhn+1∥02≤Cτ∥vhn+1∥02+Ch4×∫tn-1tn+1‍(∥ut(t)∥22+∥pt(t)∥12)dt.
The definition E^hn+1=Ehn+1+∇δψhn+1 and Assumption 1 give us
(131)A3=4μτ〈∇δψhn+1,∇rhn+1〉≤μτ∥Ehn+1∥02+Cτμ∥∇δψhn+1∥02.
If we apply (127) again, then we arrive at
(132)A3≤μτ∥Ehn+1∥02+Cτ3μ∥∇δεhn+1∥02+Cτ4μ∫tntn+1‍∥∇pt(t)∥02dt.
On the other hand, the truncation error term becomes
(133)A4=4τ〈Rn+1,vhn+1〉≤Cτ∥∇vhn+1∥02+Cτ4×∫tn-1tn+1‍∥uttt(s)∥02dt.
For the new term A5, we observe ∥δδθ(tn+1)∥0≤Cτ3∫tn-1tn+1‍∥θt(t)∥02dt; whence
(134)A5=4κμ2τ×〈g(δδθ(tn+1)+2ϑn-ϑn-1),vhn+1〉≤Cτ∥∇vhn+1∥02+Cκ2μ4τ×∥2ϑn-ϑn-1∥02+Cκ2μ4τ4×∫tn-1tn+1‍∥θt(t)∥02dt.
Invoking v0=0 and inserting the above estimates from A1 and A5 into (121) lead us to
(135)δ∥∇vhn+1∥02+δ∥∇(2vhn+1-vhn)∥02+∥∇δδvhn+1∥02+2μτ∥Ehn+1∥02≤μτ2(∥Ehn∥02+∥Gn+1∥02+∥Gn∥02+∥δEhn∥02+∥δGn∥02)+Cκ2μ4τ∥2ϑn-ϑn-1∥02+Cμτ(τ+h)2×(∥∇E^hn+1∥02+∥∇Gn+1∥02)+Cτ3μ∥∇δεhn+1∥02+Cμ(τ4+h4)×∫tn-1tn+1‍(∥ut(t)∥22+∥pt(t)∥12)dt+Cτ∥∇vhn+1∥02+Cτ4×∫tn-1tn+1(∥utt(t)∥02+∥uttt(t)∥02+∥θt(t)∥02)dt.‍
On the other hand, we choose ϕh=4τωhn+1 in (83) to obtain
(136)δ∥∇ωhn+1∥02+δ∥∇(2ωhn+1-ωhn)∥02+∥∇δδωhn+1∥02+4λμτ∥ϑhn+1∥02=A6+A7,
where
(137)A6:=4τ𝒩(uhn+1,θhn+1,ωhn+1)-4τ𝒩(u(tn+1),θ(tn+1),ωhn+1),A7∶=-2〈3ηn+1-4ηn+ηn-1,ωhn+1〉+4τ〈Qn+1,ωhn+1〉.
In order to estimate A6, we note first ∥En+1∥L3(Ω)≤Ch-1/2∥En+1∥0 and ∥ϑn+1∥12≤C((τ2+h2)/τ) which come from Lemmas 5 and 13, respectively. Then Lemma 10 and the assumption τ≤Ch yield
(138)A6=-4τ𝒩(En+1,θ(tn+1),ωhn+1)+4τ𝒩(En+1-u(tn+1),ϑn+1,ωhn+1)≤Cτ∥En+1∥0∥θ(tn+1)∥2∥ωhn+1∥1+Cτ∥En+1∥L3(Ω)∥ϑn+1∥1∥ωhn+1∥1+Cτ∥u(tn+1)∥2∥ϑn+1∥0∥ωhn+1∥1≤Cτ∥∇ωhn+1∥02+μτ∥En+1∥02+λμτ(∥ϑhn+1∥02+∥ηn+1∥02).
Since ∥δηn+1∥02≤Ch4τ∫nn+1‍∥θt(t)∥22dt, we can readily get
(139)A7≤Cτ∥∇ωhn+1∥02+C(τ4+h4)×∫n-1n+1‍(∥θt(t)∥22+∥θttt(t)∥22)dt.
Inserting the above estimates into (136) yields
(140)δ∥∇ωhn+1∥02+δ∥∇(2ωhn+1-ωhn)∥02+∥∇δδωhn+1∥02+3λμτ∥ϑhn+1∥02≤Cτ∥∇ωhn+1∥02+λμτ∥ηn+1∥02+μτ∥En+1∥02+C(τ4+h4)×∫n-1n+1‍(∥θt(t)∥22+∥θttt(t)∥22)dt.
In conjunction with the discrete Gronwall inequality, adding (135) and (140) and summing over n from 1 to N lead to (120).We now estimate the pressure error inL2(0,T;L2(Ω)). 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 τ≤Ch. If Assumptions 1–5 hold, then one has
(141)τ∑n=1N‍∥ehn+1∥02≤C(τ2+h2).Proof.
We first recall again inf-sup condition in Assumption2. Consequently, it suffices to estimate 〈en+1,∇·w〉 in terms of ∥∇w∥0. In conjunction with (10), we can rewrite (66) as
(142)〈ehn+1,∇·wh〉=12τ〈3En+1-4En+En-1,wh〉+μ〈∇E^hn+1,∇wh〉+𝒩(δδu(tn+1),u(tn+1),wh)+𝒩(2u(tn+1)-u(tn),E^n+1,wh)+𝒩(2En-En-1,u^hn+1,wh)-μ〈∇δqhn+1,wh〉-〈Rn+1,wh〉+κμ2〈g(θ(tn+1)-2θhn+θhn-1),wh〉:=∑i=18‍Ai.
We now proceed to estimate each term from A1 to A7 separately. We readily obtain
(143)A1≤Cτ(∥δEn+1∥0+∥δEn∥0)∥wh∥0≤Cτ(∥δEn+1∥0+∥δEn∥0)∥∇wh∥0,A2≤C∥∇E^hn+1∥0∥∇wh∥0.
Terms A3 and A4 can be dealt with thanks to the aid of Lemma 10 and ∥u(tn+1)∥2≤M as follows:
(144)A3≤C∥δδu(tn+1)∥0×∥u(tn+1)∥2∥wh∥1≤C∥δδu(tn+1)∥1∥∇wh∥1,A4≤C∥2u(tn+1)-u(tn)∥2×∥E^n+1∥0∥wh∥1≤C∥E^n+1∥0∥∇wh∥0.
In light of ∥u^hn+1∥1=∥E^n+1-u(tn+1)∥1≤C from Lemma 13, we can have
(145)A5≤C∥2En-En-1∥L3(Ω)×∥u^hn+1∥1∥wh∥1≤Ch(∥En+1∥0+∥En∥0)×∥∇wh∥0.
Integrate by parts and Hölder inequality yield
(146)A6≤C∥δqhn+1∥0∥∇wh∥0.
On the other hand, we have
(147)A7≤∥Rn+1∥-1∥∇w∥0.
The new term A8 can be bound by the Hölder inequality as
(148)A8=κμ2×〈g(δδθ(tn+1)-2ϑn+ϑn-1),whn+1〉≤Cκμ2(∥δδθ(tn+1)∥0+∥2ϑhn-ϑhn-1∥0+∥2ηn-ηn-1∥0)∥∇w∥0.
Inserting the estimates from A1 to A7 back into (142) and employing the discrete inf-sup condition in Assumption 2, we obtain
(149)C∥en+1∥0≤1τ(∥δEn+1∥0+∥δEn∥0)+1h(∥En+1∥0+∥En∥0)+∥δqhn+1∥0+∥∇E^hn+1∥0+∥δδu(tn+1)∥1+∥E^n+1∥0+∥Rn+1∥-1+∥δδθ(tn+1)∥0+∥2ϑhn-ϑhn-1∥0+∥2ηn-ηn-1∥0.
If we now square it, multiply it by τ, and sum over n from 1 to N, 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
(150)u=π(t2-t+1)x2(1-x)2sin(2πy),v=-2(t2-t+1)x(2x-1)(x-1)sin2(πy),p=-(t2-t-1)cos(πx)(y2+1),θ=cos(t)sin(πx)y(1-y).
Table 1 is error decay with μ=1 and τ=h. We conclude that the numerical accuracy of SGUM is optimal and consists with the result of Theorem 4.Table 1
Error decay for Algorithm1 with τ=h and μ=1.
τ
=
h
1/16
1/32
1/64
1/128
1/256
∥
E
∥
0
0.00198001
0.000651734
0.000188113
5.0822
e
-
05
1.32686
e
-
05
Order
1.603153
1.792684
1.888075
1.937437
∥
E
∥
L
∞
0.00837795
0.00282504
0.000825997
0.000225008
5.90462
e
-
05
Order
1.568326
1.774063
1.876160
1.930060
∥
E
∥
1
0.0279984
0.00950852
0.00280233
0.00077052
0.000205078
Order
1.558052
1.762594
1.862723
1.909660
∥
e
∥
0
0.0341734
0.0123334
0.00376493
0.00105695
0.000283961
Order
1.470303
1.711876
1.832716
1.896142
∥
e
∥
L
∞
0.197857
0.107171
0.0453294
0.0170497
0.00610541
Order
0.884544
1.241396
1.410701
1.481586
∥
ϑ
∥
0
0.000215482
5.3851
e
-
05
1.34612
e
-
05
3.36511
e
-
06
8.4091
e
-
07
Order
2.000522
2.000166
2.000081
2.000630
∥
ϑ
∥
L
∞
0.000197687
4.98479
e
-
05
1.25012
e
-
05
3.1314
e
-
06
7.84502
e
-
07
Order
1.987613
1.995466
1.997187
1.996959
∥
ϑ
∥
1
0.0167834
0.00839145
0.0041957
0.00209784
0.00104892
Order
1.000043
1.000009
1.000007
1.000000Example 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 f=0 and b=0. Figure 1 displays the initial and boundary conditions for velocity u and temperature θ, as already studied in [2]. Figures 2–4 are simulations at t=1 with the nondimensional parameters κ=104,λ=1, and  μ=1. Figure 2 is the result for the case τ=10-2,h=2-4 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-3,h=2-4, and τ=10-4,h=2-6, 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
Example20: initial and boundary values of Bénard convection problem.Figure 2
Example20: streamlines of velocity and isolines of temperature and pressure, at time t=1.0. The nondimensional parameters are κ=104,  λ=1,  μ=1, and the discretization parameters are τ=10-2,  h=2-4.Figure 3
Example20: streamlines of velocity and isolines of temperature and pressure, at time t=1.0. The nondimensional parameters are κ=104,λ=1, and μ=1, and the discretization parameters are τ=10-3,h=2-4.Figure 4
Example20: streamlines of velocity and isolines of temperature and pressure, at time t=1.0. The nondimensional parameters are κ=104,λ=1, and μ=1, and the discretization parameters are τ=10-4,h=2-6.Example 21 (thermal driven cavity flow).
We consider the thermal driven cavity flow in an enclosed squareΩ=[0,1]2, 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 (t=0) to steady state (t=0.2).Figure 5
Example21: initial and boundary values for thermal driven cavity flow.Example21: time sequence t=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,h=2-5. Note that umax stands for ∥u∥∞.
(a)
(b)
(c)
(d)



---

*Source: 372906-2013-05-09.xml*