# The Crank-Nicolson Extrapolation Stabilized Finite Element Method for Natural Convection Problem

Authors: Yunzhang Zhang; Yanren Hou
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

This paper studies a fully discrete Crank-Nicolson linear extrapolation stabilized
finite element method for the natural convection problem, which is
unconditionally stable and has second order temporal accuracy ofO
(
Δ
t
2
+
h
Δ
t
+
h
m
). A simple artificial viscosity stabilized of the linear system for the
approximation of the new time level connected to antidiffusion of its effects
at the old time level is used. An unconditionally stability and an a priori
error estimate are derived for the fully discrete scheme. A series of numerical
results are presented that validate our theoretical findings.


---

## Body

## 1. Introduction

Natural convection flow has many thermal engineering applications such as in double-glazed windows, solar collectors, cooling devices for electronic instruments, gas-filled cavities around nuclear reactor cores, and building insulation. Typically, fluid flow and heat transfer are governed by the partial differential equation system of momentum, mass, and energy conservation, but in the case of natural convection, the so-called Boussinesq approximation is generally used. The natural convection problem which we consider is for bounded, polyhedral domainsΩ
e
⊂
Ω in R
d (d
=
2,3) with dist
⁡
(
∂
Ω
e
,
∂
Ω
)
>
0, the simulation time t
*, and the force field γ
:
Ω
×
(
0
,
t
*
]
→
R; find the velocity u
:
Ω
×
(
0
,
t
*
]
→
R
d, the pressure p
:
Ω
×
(
0
,
t
*
]
→
R, and the temperature T
:
Ω
×
(
0
,
t
*
]
→
R satisfying [1]
(1)
u
t
-
Pr
Δ
u
+
(
u
·
∇
)
u
+
∇
p
=
Pr
Ra
ζ
T
,
ζ
=
g
|
g
|
,
u
=
0
on
⁡
∂
Ω
e
,
u
≡
0
in
Ω
-
Ω
e
=
Ω
s
,
u
|
t
=
0
=
u
0
,
∇
·
u
(
x
,
t
)
=
0
in
Ω
e
,
T
t
-
∇
·
(
k
∇
T
)
+
(
u
·
∇
)
T
=
γ
in
Ω
,
T
=
0
on
⁡
Γ
T
,
∂
T
∂
n
=
0
on
⁡
Γ
B
,
T
|
t
=
0
=
T
0
,
in
Ω
,
where ζ is a unit vector in the direction of gravity, n is the outward unit normal to Ω, and Γ
T
=
∂
Ω
∖
Γ
B where Γ
B is a regular open subset of ∂
Ω, Pr is Prandtl number, Ra is Rayleigh number, and k
>
0 is thermal conductivity parameter. Moreover, k
=
k
e in Ω
e and k
=
k
s in Ω
s, where k
e and k
s are positive constants. A global-in-time existence result for a more general natural convection problem (Navier-Stokes/Fourier model) is given in [2].Many authors have worked hard to study for a great variety of efficient numerical schemes for the natural convection problem [3–17] and relevant research [18, 19]. We mention only a few papers here. [3, 4] are the early papers by using mixed finite element (FE) method. Çıbık and Kaya [5] have formulated a projection-based stabilization FE technique for solving the steady-state natural convection problems. The global stabilizations are added for both velocity and temperature variables and these effects are subtracted from the large scales. Galvin et al. [7] consider the problem of poor mass conservation in mixed FE algorithms for flow problems with large rotation-free forcing in the momentum equation. Zhang et al. [8] have presented a subgrid stabilized defect-correction method for steady-state natural convection problem. Shi and Ren [11] have proposed a least squares Galerkin-Petrov nonconforming mixed FE method for stationary conduction-convection problems. Luo et al. [12] have given an optimizing reduced Petrov-Galerkin least squares mixed FE for the nonstationary conduction-convection problem. Boland and Layton [1] have derived stability properties and error estimates for the mixed FE spatial discretization case when used to approximate heat flow in a fluid enclosed by a solid medium. Benítez and Bermúdez have presented a second order Lagrange-Galerkin method for natural convection problems in [17]. In [20, 21], a stability analysis of thermal natural convection in superposed fluid and porous layers is carried out.Our goal in this paper is to solve time-dependent natural convection problem efficiently and accurately. Usually fully implicit schemes are (almost) unconditionally stable, but one has to solve a system of nonlinear equations at each time step. Although an explicit scheme is much easier in computation, it suffers a restricted time step size from the stability requirement. A popular approach is based on an implicit scheme for the linear term and a semi-implicit scheme or an explicit scheme for the nonlinear term. There are numerous works on the Crank-Nicolson and relevant high order scheme for the Navier-Stokes (NS) equations [22–30]. The Crank-Nicolson linear extrapolation (CNLE) scheme for NS equations was first studied by Baker in [23]. The second and third order CNLE methods are introduced and analysed in [24]. A stabilized extrapolated trapezoidal FE method is given in [25] for the NS equations. A variational multiscale method based on the CNLE scheme for the NS equations is proposed in [26]. He et al. [27–29] have studied the NS equations based on the Crank-Nicolson extrapolation (Crank-Nicolson/Adams-Bashforth, or two level methods) schemes. In [31], we have studied fully implicit Crank-Nicolson scheme for natural convection problem.We consider herein a simple, second order accurate, and unconditionally stable fully discrete Crank-Nicolson linear extrapolation stabilized (CNSLE) FE method for natural convection problem which requires the solution of one linear system per time step. Suppressing the spatial discretization, the method is(2)
u
n
+
1
-
u
n
Δ
t
-
Pr
Δ
(
u
n
+
1
+
u
n
2
)
+
(
U
n
+
1
/
2
·
∇
)
(
u
n
+
1
+
u
n
2
)
+
(
p
n
+
1
+
p
n
2
)
-
μ
h
Δ
u
n
+
1
=
Pr
Ra
ζ
(
T
n
+
1
+
T
n
2
)
-
μ
h
Δ
u
n
,
∇
·
u
n
+
1
=
0
,
T
n
+
1
-
T
n
Δ
t
-
∇
·
(
k
∇
(
T
n
+
1
+
T
n
2
)
)
+
(
U
n
+
1
/
2
·
∇
)
(
T
n
+
1
+
T
n
2
)
-
μ
h
Δ
T
n
+
1
=
γ
n
+
1
/
2
-
μ
h
Δ
T
n
,
where the time step Δ
t
>
0, the constant μ
=
O
(
1
), and U
n
+
1
/
2
=
(
3
/
2
)
u
n
-
(
1
/
2
)
u
n
-
1 is the linear extrapolation of the velocity to t
n
+
1
/
2 from previous time levels. It is a three time levels scheme. Artificial viscosity stabilizations are introduced into the linear systems for u
n
+
1 and T
n
+
1 by adding -
μ
h
Δ
u
n
+
1 and -
μ
h
Δ
T
n
+
1 to the left-hand sides (LHS) and correcting them by -
μ
h
Δ
u
n and -
μ
h
Δ
T
n on the right-hand sides (RHS), respectively. To the best of the authors’ knowledge, there are no papers dealing with the error analysis of the fully discrete CNLE FE method for natural convection problem.The paper is organized as follows. Section2 collects some preliminaries for the analysis that follows. In Section 3, we give the fully discrete CNSLE FE method and prove it is unconditionally stable in Theorem 5. In Section 4, error estimates for velocity and temperature are derived in Theorem 6. Numerical experiments are described in Section 5. Conclusions follow.


## 2. Preliminaries

### 2.1. The Variational Formulation of Natural Convection Problem

Let(
·
,
·
) and ∥
·
∥ denote the L
2
(
Ω
) inner product and norm, respectively. Define the velocity space X, the pressure space M, the temperature space W, and the divergence-free space V as follows:
(3)
X
∶
=
H
0
1
(
Ω
)
d
=
{
v
∈
H
1
(
Ω
)
d
:
v
=
0
on
⁡
Γ
}
,
M
∶
=
L
0
2
(
Ω
)
=
{
p
∈
L
2
(
Ω
)
,
(
p
,
1
)
Ω
=
0
}
,
W
∶
=
{
S
∈
H
1
(
Ω
)
:
S
=
0
on
⁡
Γ
}
,
V
∶
=
H
0
,
div
⁡
1
(
Ω
)
=
{
v
∈
X
:
∇
·
v
=
0
in
Ω
}
,
where H
j
(
Ω
) denotes the standard Sobolev space [32] with norm ∥
·
∥
j. All other norms will be clearly labeled with subscripts.The weak formulation of problem (1) reads as follows: find (
u
,
p
,
T
)
∈
(
X
,
M
,
W
) for all t
∈
(
0
,
t
*
], for all (
v
,
q
,
S
)
∈
(
X
,
M
,
W
), such that
(4)
(
u
t
,
v
)
+
Pr
(
∇
u
,
∇
v
)
+
c
(
u
,
u
,
v
)
-
(
∇
·
v
,
p
)
=
Pr
Ra
(
ζ
T
,
v
)
,
(
∇
·
u
,
q
)
=
0
,
(
T
t
,
S
)
+
k
(
∇
T
,
∇
S
)
+
c
-
(
u
,
T
,
S
)
=
(
γ
,
S
)
.
Here, the skew-symmetric trilinear forms [1, 3, 5]
(5)
c
(
u
,
v
,
w
)
=
1
2
(
(
u
·
∇
)
v
·
w
)
-
1
2
(
(
u
·
∇
)
w
·
v
)
,
vvvvvvvvvvvvvvvvv
∀
u
,
v
,
w
∈
X
,
c
-
(
u
,
T
,
S
)
=
1
2
(
(
u
·
∇
)
T
S
)
-
1
2
(
(
u
·
∇
)
S
T
)
,
vvvvvvvv
∀
u
∈
X
;
T
,
S
∈
W
,
which satisfy the following lemma.Lemma 1 (see [25, 33]).
LetΩ
⊂
R
d (d
=
2,3), for all u
,
v
,
w
∈
X; T
,
S
∈
W. 									(6)
c
(
u
,
v
,
w
)
=
-
c
(
u
,
w
,
v
)
,
|
c
(
u
,
v
,
w
)
|
⩽
C
∥
∇
u
∥
∥
∇
v
∥
∥
∇
w
∥
,
|
c
(
u
,
v
,
w
)
|
⩽
C
∥
u
∥
1
/
2
∥
∇
u
∥
1
/
2
∥
∇
v
∥
∥
∇
w
∥
,
c
-
(
u
,
T
,
S
)
=
-
c
-
(
u
,
S
,
T
)
,
c
-
(
u
,
T
,
S
)
⩽
C
∥
∇
u
∥
∥
∇
T
∥
∥
∇
S
∥
,
c
-
(
u
,
T
,
S
)
⩽
C
∥
u
∥
1
/
2
∥
∇
u
∥
1
/
2
∥
∇
T
∥
∥
∇
S
∥
,
if, in addition, v
,
∇
v
∈
L
∞
(
Ω
),
(7)
|
c
(
u
,
v
,
w
)
|
⩽
C
(
∥
v
∥
L
∞
(
Ω
)
+
∥
∇
v
∥
L
∞
(
Ω
)
)
∥
∇
u
∥
∥
∇
w
∥
,
|
c
(
u
,
v
,
w
)
|
⩽
C
(
∥
u
∥
∥
∇
v
∥
L
∞
(
Ω
)
+
∥
∇
u
∥
∥
v
∥
L
∞
(
Ω
)
)
∥
w
∥
,
and if, in addition, T
,
∇
T
∈
L
∞
(
Ω
),
(8)
|
c
-
(
u
,
T
,
S
)
|
⩽
C
(
∥
T
∥
L
∞
(
Ω
)
+
∥
∇
T
∥
L
∞
(
Ω
)
)
∥
∇
u
∥
∥
∇
S
∥
,
|
c
-
(
u
,
T
,
S
)
|
⩽
C
(
∥
u
∥
∥
∇
T
∥
L
∞
(
Ω
)
+
∥
∇
u
∥
∥
T
∥
L
∞
(
Ω
)
)
∥
S
∥
.We will use the Poincare inequality: for allω
∈
X, or W, ∥
ω
∥
⩽
C
P
∥
∇
ω
∥.


### 2.2. Finite Element Approximation

AssumeΩ
h
=
{
K
} to be a quasiuniform mesh of Ω with mesh size 0
<
h
<
1; let (
X
h
,
M
h
,
W
h
)
⊂
(
X
,
M
,
W
) be a pair of conforming velocity-pressure-temperature finite element spaces which contain piecewise continuous polynomials of degree m, m
-
1, and m, respectively, and satisfy the usual inf-sup condition and the following approximation properties [33]:
(9)
inf
⁡
v
h
∈
X
h
{
∥
u
-
v
h
∥
+
h
∥
∇
(
u
-
v
h
)
∥
}
≤
C
h
m
+
1
∥
u
∥
m
+
1
,
vvvvvvvvivvvvvvvvvvvvvv
∀
u
∈
H
m
+
1
(
Ω
)
∩
X
,
inf
⁡
q
h
∈
M
h
∥
p
-
q
h
∥
≤
C
h
m
∥
p
∥
m
,
∀
p
∈
H
m
(
Ω
)
∩
M
,
inf
⁡
S
h
∈
W
h
{
∥
T
-
S
h
∥
+
h
∥
∇
(
T
-
S
h
)
∥
}
≤
C
h
m
+
1
∥
T
∥
m
+
1
,
vvvvvvvvvvvvivvvvvvvvvv
∀
T
∈
H
m
+
1
(
Ω
)
∩
W
.
The subspace V
h of X
h is given by
(10)
V
h
=
{
v
h
∈
X
h
:
(
q
h
,
∇
·
v
h
)
=
0
,
∀
q
h
∈
M
h
}
.


### 2.3. The Modified Stokes Projection

The following modified Stokes projection operators are similar than the one in [5]. For the reader’s convenience, we only present the definition and the error results of the projection operators. We can easily derive the result.Definition 2 (see [5]).
The modified Stokes projection operator for velocityu and pressure p is P
S
:
(
X
,
M
)
→
(
X
h
,
M
h
), P
S
(
u
,
p
)
=
(
u
~
,
p
~
), such that
(11)
Pr
(
∇
(
u
-
u
~
)
,
∇
v
h
)
-
(
p
-
p
~
,
∇
·
v
h
)
=
0
,
(
∇
·
(
u
-
u
~
)
,
q
h
)
=
0
,
for all (
v
h
,
q
h
)
∈
(
X
h
,
M
h
). In spaces (V
h
,
M
h), Definition 2 reduces to the following: given (
u
,
p
)
∈
(
X
,
M
), find u
~
∈
V
h such that
(12)
Pr
(
∇
(
u
-
u
~
)
,
∇
v
h
)
-
(
p
-
q
h
,
∇
·
v
h
)
=
0
,
for all (
v
h
,
q
h
)
∈
(
V
h
,
M
h
). The modified Stokes projection operator for temperature T is P
T
:
W
→
W
h, P
T
(
T
)
=
T
~, such that
(13)
(
∇
(
T
-
T
~
)
,
∇
S
h
)
=
0
,
∀
S
h
∈
W
h
.Lemma 3 (error estimates of the modified Stokes projection).
Suppose the inf-sup condition holds; then(
u
~
,
T
~
) exists uniquely and satisfies
(14)
P
r
∥
∇
(
u
-
u
~
)
∥
2
⩽
C
(
P
r
inf
⁡
v
h
∈
X
h
∥
∇
(
u
-
v
h
)
∥
2
+
P
r
-
1
inf
⁡
q
h
∈
M
h
∥
p
-
q
h
∥
2
)
,
∥
∇
(
T
-
T
~
)
∥
⩽
C
inf
⁡
S
h
∈
W
h
∥
∇
(
T
-
S
h
)
∥
,
where C is a constant independent of h and m.


## 2.1. The Variational Formulation of Natural Convection Problem

Let(
·
,
·
) and ∥
·
∥ denote the L
2
(
Ω
) inner product and norm, respectively. Define the velocity space X, the pressure space M, the temperature space W, and the divergence-free space V as follows:
(3)
X
∶
=
H
0
1
(
Ω
)
d
=
{
v
∈
H
1
(
Ω
)
d
:
v
=
0
on
⁡
Γ
}
,
M
∶
=
L
0
2
(
Ω
)
=
{
p
∈
L
2
(
Ω
)
,
(
p
,
1
)
Ω
=
0
}
,
W
∶
=
{
S
∈
H
1
(
Ω
)
:
S
=
0
on
⁡
Γ
}
,
V
∶
=
H
0
,
div
⁡
1
(
Ω
)
=
{
v
∈
X
:
∇
·
v
=
0
in
Ω
}
,
where H
j
(
Ω
) denotes the standard Sobolev space [32] with norm ∥
·
∥
j. All other norms will be clearly labeled with subscripts.The weak formulation of problem (1) reads as follows: find (
u
,
p
,
T
)
∈
(
X
,
M
,
W
) for all t
∈
(
0
,
t
*
], for all (
v
,
q
,
S
)
∈
(
X
,
M
,
W
), such that
(4)
(
u
t
,
v
)
+
Pr
(
∇
u
,
∇
v
)
+
c
(
u
,
u
,
v
)
-
(
∇
·
v
,
p
)
=
Pr
Ra
(
ζ
T
,
v
)
,
(
∇
·
u
,
q
)
=
0
,
(
T
t
,
S
)
+
k
(
∇
T
,
∇
S
)
+
c
-
(
u
,
T
,
S
)
=
(
γ
,
S
)
.
Here, the skew-symmetric trilinear forms [1, 3, 5]
(5)
c
(
u
,
v
,
w
)
=
1
2
(
(
u
·
∇
)
v
·
w
)
-
1
2
(
(
u
·
∇
)
w
·
v
)
,
vvvvvvvvvvvvvvvvv
∀
u
,
v
,
w
∈
X
,
c
-
(
u
,
T
,
S
)
=
1
2
(
(
u
·
∇
)
T
S
)
-
1
2
(
(
u
·
∇
)
S
T
)
,
vvvvvvvv
∀
u
∈
X
;
T
,
S
∈
W
,
which satisfy the following lemma.Lemma 1 (see [25, 33]).
LetΩ
⊂
R
d (d
=
2,3), for all u
,
v
,
w
∈
X; T
,
S
∈
W. 									(6)
c
(
u
,
v
,
w
)
=
-
c
(
u
,
w
,
v
)
,
|
c
(
u
,
v
,
w
)
|
⩽
C
∥
∇
u
∥
∥
∇
v
∥
∥
∇
w
∥
,
|
c
(
u
,
v
,
w
)
|
⩽
C
∥
u
∥
1
/
2
∥
∇
u
∥
1
/
2
∥
∇
v
∥
∥
∇
w
∥
,
c
-
(
u
,
T
,
S
)
=
-
c
-
(
u
,
S
,
T
)
,
c
-
(
u
,
T
,
S
)
⩽
C
∥
∇
u
∥
∥
∇
T
∥
∥
∇
S
∥
,
c
-
(
u
,
T
,
S
)
⩽
C
∥
u
∥
1
/
2
∥
∇
u
∥
1
/
2
∥
∇
T
∥
∥
∇
S
∥
,
if, in addition, v
,
∇
v
∈
L
∞
(
Ω
),
(7)
|
c
(
u
,
v
,
w
)
|
⩽
C
(
∥
v
∥
L
∞
(
Ω
)
+
∥
∇
v
∥
L
∞
(
Ω
)
)
∥
∇
u
∥
∥
∇
w
∥
,
|
c
(
u
,
v
,
w
)
|
⩽
C
(
∥
u
∥
∥
∇
v
∥
L
∞
(
Ω
)
+
∥
∇
u
∥
∥
v
∥
L
∞
(
Ω
)
)
∥
w
∥
,
and if, in addition, T
,
∇
T
∈
L
∞
(
Ω
),
(8)
|
c
-
(
u
,
T
,
S
)
|
⩽
C
(
∥
T
∥
L
∞
(
Ω
)
+
∥
∇
T
∥
L
∞
(
Ω
)
)
∥
∇
u
∥
∥
∇
S
∥
,
|
c
-
(
u
,
T
,
S
)
|
⩽
C
(
∥
u
∥
∥
∇
T
∥
L
∞
(
Ω
)
+
∥
∇
u
∥
∥
T
∥
L
∞
(
Ω
)
)
∥
S
∥
.We will use the Poincare inequality: for allω
∈
X, or W, ∥
ω
∥
⩽
C
P
∥
∇
ω
∥.


## 2.2. Finite Element Approximation

AssumeΩ
h
=
{
K
} to be a quasiuniform mesh of Ω with mesh size 0
<
h
<
1; let (
X
h
,
M
h
,
W
h
)
⊂
(
X
,
M
,
W
) be a pair of conforming velocity-pressure-temperature finite element spaces which contain piecewise continuous polynomials of degree m, m
-
1, and m, respectively, and satisfy the usual inf-sup condition and the following approximation properties [33]:
(9)
inf
⁡
v
h
∈
X
h
{
∥
u
-
v
h
∥
+
h
∥
∇
(
u
-
v
h
)
∥
}
≤
C
h
m
+
1
∥
u
∥
m
+
1
,
vvvvvvvvivvvvvvvvvvvvvv
∀
u
∈
H
m
+
1
(
Ω
)
∩
X
,
inf
⁡
q
h
∈
M
h
∥
p
-
q
h
∥
≤
C
h
m
∥
p
∥
m
,
∀
p
∈
H
m
(
Ω
)
∩
M
,
inf
⁡
S
h
∈
W
h
{
∥
T
-
S
h
∥
+
h
∥
∇
(
T
-
S
h
)
∥
}
≤
C
h
m
+
1
∥
T
∥
m
+
1
,
vvvvvvvvvvvvivvvvvvvvvv
∀
T
∈
H
m
+
1
(
Ω
)
∩
W
.
The subspace V
h of X
h is given by
(10)
V
h
=
{
v
h
∈
X
h
:
(
q
h
,
∇
·
v
h
)
=
0
,
∀
q
h
∈
M
h
}
.


## 2.3. The Modified Stokes Projection

The following modified Stokes projection operators are similar than the one in [5]. For the reader’s convenience, we only present the definition and the error results of the projection operators. We can easily derive the result.Definition 2 (see [5]).
The modified Stokes projection operator for velocityu and pressure p is P
S
:
(
X
,
M
)
→
(
X
h
,
M
h
), P
S
(
u
,
p
)
=
(
u
~
,
p
~
), such that
(11)
Pr
(
∇
(
u
-
u
~
)
,
∇
v
h
)
-
(
p
-
p
~
,
∇
·
v
h
)
=
0
,
(
∇
·
(
u
-
u
~
)
,
q
h
)
=
0
,
for all (
v
h
,
q
h
)
∈
(
X
h
,
M
h
). In spaces (V
h
,
M
h), Definition 2 reduces to the following: given (
u
,
p
)
∈
(
X
,
M
), find u
~
∈
V
h such that
(12)
Pr
(
∇
(
u
-
u
~
)
,
∇
v
h
)
-
(
p
-
q
h
,
∇
·
v
h
)
=
0
,
for all (
v
h
,
q
h
)
∈
(
V
h
,
M
h
). The modified Stokes projection operator for temperature T is P
T
:
W
→
W
h, P
T
(
T
)
=
T
~, such that
(13)
(
∇
(
T
-
T
~
)
,
∇
S
h
)
=
0
,
∀
S
h
∈
W
h
.Lemma 3 (error estimates of the modified Stokes projection).
Suppose the inf-sup condition holds; then(
u
~
,
T
~
) exists uniquely and satisfies
(14)
P
r
∥
∇
(
u
-
u
~
)
∥
2
⩽
C
(
P
r
inf
⁡
v
h
∈
X
h
∥
∇
(
u
-
v
h
)
∥
2
+
P
r
-
1
inf
⁡
q
h
∈
M
h
∥
p
-
q
h
∥
2
)
,
∥
∇
(
T
-
T
~
)
∥
⩽
C
inf
⁡
S
h
∈
W
h
∥
∇
(
T
-
S
h
)
∥
,
where C is a constant independent of h and m.


## 3. Numerical Scheme and Its Stability

### 3.1. Numerical Scheme

Lett
n
=
n
Δ
t, n
=
0,1
,
2
,
…
,
N, and t
*
=
N
Δ
t. Define the linear extrapolation of the convecting velocity u
h to t
n
+
1
/
2
∶
=
(
t
n
+
t
n
+
1
)
/
2 by
(15)
E
[
u
h
n
,
u
h
n
-
1
]
=
3
2
u
h
n
-
1
2
u
h
n
-
1
,
where u
h
j is a known approximation to u
(
x
,
t
j
). The fully discrete CNSLE FE method of (1) is presented as follows.Algorithm 4.
Consider the following steps.
Step
1. Let (
u
h
0
,
T
h
0
) be the modified Stokes projections of (
u
0
,
T
0
) into spaces (
V
h
,
W
h
); then at the first time level, find (
u
h
1
,
p
h
1
,
T
h
1
)
∈
(
X
h
,
M
h
,
W
h
), for all (
v
h
,
q
h
,
S
h
)
∈
(
X
h
,
M
h
,
W
h
), such that
(16)
(
u
h
1
-
u
h
0
Δ
t
,
v
h
)
+
Pr
⁡
(
∇
(
u
h
1
+
u
h
0
2
)
,
∇
v
h
)
+
μ
h
(
∇
u
h
1
,
∇
v
h
)
-
(
p
h
1
+
p
h
0
2
,
∇
·
v
h
)
+
c
(
u
h
0
,
u
h
1
+
u
h
0
2
,
v
h
)
=
Pr
Ra
(
ζ
(
T
h
1
+
T
h
0
2
)
,
v
h
)
+
μ
h
(
∇
u
h
0
,
∇
v
h
)
,
(
∇
·
u
h
1
,
q
h
)
=
0
,
(17)
(
T
h
1
-
T
h
0
Δ
t
,
S
h
)
+
k
(
∇
(
T
h
1
+
T
h
0
2
)
,
∇
S
h
)
+
μ
h
(
∇
T
h
1
,
∇
S
h
)
+
c
-
(
u
h
0
,
T
h
1
+
T
h
0
2
,
S
h
)
=
(
γ
(
t
1
/
2
)
,
S
h
)
+
μ
h
(
∇
T
h
0
,
∇
S
h
)
.
Step
2. For n
⩾
1, given (
u
h
n
,
p
h
n
,
T
h
n
)
∈
(
X
h
,
M
h
,
W
h
), find (
u
h
n
+
1
,
p
h
n
+
1
,
T
h
n
+
1
)
∈
(
X
h
,
M
h
,
W
h
), for all (
v
h
,
q
h
,
S
h
)
∈
(
X
h
,
M
h
,
W
h
), such that
(18)
(
u
h
n
+
1
-
u
h
n
Δ
t
,
v
h
)
+
Pr
(
∇
(
u
h
n
+
1
+
u
h
n
2
)
,
∇
v
h
)
+
μ
h
(
∇
u
h
n
+
1
,
∇
v
h
)
+
c
(
E
[
u
h
n
,
u
h
n
-
1
]
,
u
h
n
+
1
+
u
h
n
2
,
v
h
)
-
(
p
h
n
+
1
+
p
h
n
2
,
∇
·
v
h
)
=
Pr
Ra
(
ζ
(
T
h
n
+
1
+
T
h
n
2
)
,
v
h
)
+
μ
h
(
∇
u
h
n
,
∇
v
h
)
,
(
∇
·
u
h
n
+
1
,
q
h
)
=
0
,
(19)
(
T
h
n
+
1
-
T
h
n
Δ
t
,
S
h
)
+
k
(
∇
(
T
h
n
+
1
+
T
h
n
2
)
,
∇
S
h
)
+
μ
h
(
∇
T
h
n
+
1
,
∇
S
h
)
+
c
-
(
E
[
u
h
n
,
u
h
n
-
1
]
,
T
h
n
+
1
+
T
h
n
2
,
S
h
)
=
(
γ
(
t
n
+
1
/
2
)
,
S
h
)
+
μ
h
(
∇
T
h
n
,
∇
S
h
)
.We find that CNLES FE method requires the solution of only one linear problem at each time step; thus it needs less time than fully implicit Crank-Nicolson scheme. Denotek
min
⁡
=
min
⁡
(
k
e
,
k
s
), and k
max
⁡
=
max
⁡
(
k
e
,
k
s
).


### 3.2. Stability of the Method

Theorem 5.
Suppose0
<
k
min
⁡
≤
k
max
⁡
<
∞, and γ
∈
L
2
(
0
,
T
;
H
-
1
(
Ω
)
). Algorithm 4 is unconditionally stable in the following sense, for any h, Δ
t
>
0, μ
⩾
0, and 0
⩽
l
⩽
N
-
1:(20)
∥
T
h
l
+
1
∥
2
+
k
min
⁡
Δ
t
∑
n
=
0
l
∥
∇
(
T
h
n
+
1
+
T
h
n
2
)
∥
2
+
μ
h
Δ
t
∥
∇
T
h
l
+
1
∥
2
⩽
∥
T
h
0
∥
2
+
μ
h
Δ
t
∥
∇
T
h
0
∥
2
+
k
min
⁡
-
1
Δ
t
∑
n
=
0
l
∥
γ
(
t
n
+
1
/
2
)
∥
-
1
2
,
(21)
∥
u
h
l
+
1
∥
2
+
Pr
Δ
t
∑
n
=
0
l
∥
∇
(
u
h
n
+
1
+
u
h
n
2
)
∥
2
+
μ
h
Δ
t
∥
∇
u
h
l
+
1
∥
2
⩽
∥
u
h
0
∥
2
+
μ
h
Δ
t
∥
∇
u
h
0
∥
2
+
Pr
Ra
2
k
min
⁡
-
1
C
P
[
∑
n
=
0
l
k
min
⁡
-
1
∥
γ
(
t
n
+
(
1
/
2
)
)
∥
-
1
2
∥
T
h
0
∥
2
+
μ
h
Δ
t
∥
∇
T
h
0
∥
2
0000000000000000000
+
Δ
t
∑
n
=
0
l
k
min
⁡
-
1
∥
γ
(
t
n
+
1
/
2
)
∥
-
1
2
]
.Proof.
We first derive the stability of temperatureT and then of velocity u. Choosing S
h
=
(
T
h
1
+
T
h
0
)
/
2
∈
W
h in (17) yields
(22)
(
T
h
1
-
T
h
0
Δ
t
,
T
h
1
+
T
h
0
2
)
+
k
∥
∇
(
T
h
1
+
T
h
0
2
)
∥
2
+
μ
h
Δ
t
(
∇
(
T
h
1
-
T
h
0
Δ
t
)
,
∇
(
T
h
1
+
T
h
0
2
)
)
=
(
γ
(
t
1
/
2
)
,
T
h
1
+
T
h
0
2
)
.
Applying the Cauchy-Schwarz and Young’s inequalities gives
(23)
∥
T
h
1
∥
2
-
∥
T
h
0
∥
2
2
Δ
t
+
k
min
⁡
∥
∇
(
T
h
1
+
T
h
0
2
)
∥
2
+
μ
h
Δ
t
∥
∇
T
h
1
∥
2
-
∥
T
h
0
∥
2
2
Δ
t
⩽
1
2
k
min
⁡
-
1
∥
γ
(
t
1
/
2
)
∥
-
1
2
+
1
2
k
min
⁡
∥
∇
(
T
h
1
+
T
h
0
2
)
∥
2
.
Thus, we get the bound of temperature on the first time level:
(24)
∥
T
h
1
∥
2
+
Δ
t
k
min
⁡
∥
∇
(
T
h
1
+
T
h
0
2
)
∥
2
+
μ
h
Δ
t
∥
∇
T
h
1
∥
2
⩽
∥
T
h
0
∥
2
+
μ
h
Δ
t
∥
∇
T
h
0
∥
2
+
k
min
⁡
-
1
Δ
t
∥
γ
(
t
1
/
2
)
∥
-
1
2
.
For n
⩾
1, taking S
h
=
(
T
h
n
+
1
+
T
h
n
)
/
2
∈
W
h in (19) and using the Cauchy-Schwarz and Young inequalities lead to
(25)
(
∥
T
h
n
+
1
∥
2
-
∥
T
h
n
∥
2
)
+
k
min
⁡
Δ
t
∥
∇
(
T
h
n
+
1
+
T
h
n
2
)
∥
2
+
μ
h
Δ
t
∥
∇
T
h
n
+
1
∥
2
⩽
μ
h
Δ
t
∥
∇
T
h
n
∥
2
+
k
min
⁡
-
1
Δ
t
∥
γ
(
t
n
+
1
/
2
)
∥
-
1
2
.
Summing the above equation over n from 1 to l yields
(26)
∥
T
h
l
+
1
∥
2
+
k
min
⁡
Δ
t
∑
n
=
1
l
∥
∇
(
T
h
n
+
1
+
T
h
n
2
)
∥
2
+
μ
h
Δ
t
∥
∇
T
h
l
+
1
∥
2
⩽
∥
T
h
1
∥
2
+
μ
h
Δ
t
∥
∇
T
h
1
∥
2
+
k
min
⁡
-
1
Δ
t
∑
n
=
1
l
∥
γ
(
t
n
+
1
/
2
)
∥
-
1
2
.
Substituting the bound ∥
T
h
1
∥
2
+
Δ
t
k
min
⁡
∥
∇
(
(
T
h
1
+
T
h
0
)
/
2
)
∥
2
+
μ
h
Δ
t
∥
∇
T
h
1
∥
2 of (24) into the above relation, we obtain the result of (20).
Now we derive the stability of velocityu. Taking v
h
=
(
u
h
1
+
u
h
0
)
/
2
∈
V
h in (16) and applying the Cauchy-Schwarz and Young inequalities, we have
(27)
∥
u
h
1
∥
2
-
∥
u
h
0
∥
2
2
Δ
t
+
Pr
∥
∇
(
u
h
1
+
u
h
0
2
)
∥
2
+
μ
h
Δ
t
∥
∇
u
h
1
∥
2
-
∥
∇
u
h
0
∥
2
2
Δ
t
⩽
Pr
2
∥
∇
(
u
h
1
+
u
h
0
2
)
∥
2
+
Pr
R
a
2
2
∥
T
h
1
+
T
h
0
2
∥
-
1
2
.
Then
(28)
∥
u
h
1
∥
2
+
Pr
Δ
t
∥
∇
(
u
h
1
+
u
h
0
2
)
∥
2
+
μ
h
Δ
t
∥
∇
u
h
1
∥
2
⩽
∥
u
h
0
∥
2
+
μ
h
Δ
t
∥
∇
u
h
0
∥
2
+
Pr
R
a
2
Δ
t
∥
T
h
1
+
T
h
0
2
∥
-
1
2
.
For n
⩾
1, taking v
h
=
(
u
h
n
+
1
+
u
h
n
)
/
2
∈
V
h in (18) and making use of the Cauchy-Schwarz inequality and the Young’s inequality lead to
(29)
∥
u
h
n
+
1
∥
2
-
∥
u
h
n
∥
2
2
Δ
t
+
Pr
∥
∇
(
u
h
n
+
1
+
u
h
n
2
)
∥
2
+
μ
h
Δ
t
∥
∇
u
h
n
+
1
∥
2
-
∥
∇
u
h
n
∥
2
2
Δ
t
⩽
Pr
2
∥
∇
(
u
h
n
+
1
+
u
h
n
2
)
∥
2
+
Pr
R
a
2
2
∥
T
h
n
+
1
+
T
h
n
2
∥
-
1
2
.
Multiplying the above relation by 2
Δ
t and summing over n from 1 to l give
(30)
∥
u
h
l
+
1
∥
2
+
Pr
Δ
t
∑
n
=
1
l
∥
∇
(
u
h
n
+
1
+
u
h
n
2
)
∥
2
+
μ
h
Δ
t
∥
∇
u
h
l
+
1
∥
2
⩽
∥
u
h
1
∥
2
+
μ
h
Δ
t
∥
∇
u
h
1
∥
2
+
Δ
t
Pr
Ra
2
∑
n
=
1
l
∥
T
h
n
+
1
+
T
h
n
2
∥
-
1
2
.
Inserting the bound ∥
u
h
1
∥
2
+
Pr
Δ
t
∥
∇
(
(
u
h
1
+
u
h
0
)
/
2
)
∥
2
+
μ
h
Δ
t
∥
∇
u
h
1
∥
2 of (28) into the above relation and using the result (20), we get the desired result (21).


## 3.1. Numerical Scheme

Lett
n
=
n
Δ
t, n
=
0,1
,
2
,
…
,
N, and t
*
=
N
Δ
t. Define the linear extrapolation of the convecting velocity u
h to t
n
+
1
/
2
∶
=
(
t
n
+
t
n
+
1
)
/
2 by
(15)
E
[
u
h
n
,
u
h
n
-
1
]
=
3
2
u
h
n
-
1
2
u
h
n
-
1
,
where u
h
j is a known approximation to u
(
x
,
t
j
). The fully discrete CNSLE FE method of (1) is presented as follows.Algorithm 4.
Consider the following steps.
Step
1. Let (
u
h
0
,
T
h
0
) be the modified Stokes projections of (
u
0
,
T
0
) into spaces (
V
h
,
W
h
); then at the first time level, find (
u
h
1
,
p
h
1
,
T
h
1
)
∈
(
X
h
,
M
h
,
W
h
), for all (
v
h
,
q
h
,
S
h
)
∈
(
X
h
,
M
h
,
W
h
), such that
(16)
(
u
h
1
-
u
h
0
Δ
t
,
v
h
)
+
Pr
⁡
(
∇
(
u
h
1
+
u
h
0
2
)
,
∇
v
h
)
+
μ
h
(
∇
u
h
1
,
∇
v
h
)
-
(
p
h
1
+
p
h
0
2
,
∇
·
v
h
)
+
c
(
u
h
0
,
u
h
1
+
u
h
0
2
,
v
h
)
=
Pr
Ra
(
ζ
(
T
h
1
+
T
h
0
2
)
,
v
h
)
+
μ
h
(
∇
u
h
0
,
∇
v
h
)
,
(
∇
·
u
h
1
,
q
h
)
=
0
,
(17)
(
T
h
1
-
T
h
0
Δ
t
,
S
h
)
+
k
(
∇
(
T
h
1
+
T
h
0
2
)
,
∇
S
h
)
+
μ
h
(
∇
T
h
1
,
∇
S
h
)
+
c
-
(
u
h
0
,
T
h
1
+
T
h
0
2
,
S
h
)
=
(
γ
(
t
1
/
2
)
,
S
h
)
+
μ
h
(
∇
T
h
0
,
∇
S
h
)
.
Step
2. For n
⩾
1, given (
u
h
n
,
p
h
n
,
T
h
n
)
∈
(
X
h
,
M
h
,
W
h
), find (
u
h
n
+
1
,
p
h
n
+
1
,
T
h
n
+
1
)
∈
(
X
h
,
M
h
,
W
h
), for all (
v
h
,
q
h
,
S
h
)
∈
(
X
h
,
M
h
,
W
h
), such that
(18)
(
u
h
n
+
1
-
u
h
n
Δ
t
,
v
h
)
+
Pr
(
∇
(
u
h
n
+
1
+
u
h
n
2
)
,
∇
v
h
)
+
μ
h
(
∇
u
h
n
+
1
,
∇
v
h
)
+
c
(
E
[
u
h
n
,
u
h
n
-
1
]
,
u
h
n
+
1
+
u
h
n
2
,
v
h
)
-
(
p
h
n
+
1
+
p
h
n
2
,
∇
·
v
h
)
=
Pr
Ra
(
ζ
(
T
h
n
+
1
+
T
h
n
2
)
,
v
h
)
+
μ
h
(
∇
u
h
n
,
∇
v
h
)
,
(
∇
·
u
h
n
+
1
,
q
h
)
=
0
,
(19)
(
T
h
n
+
1
-
T
h
n
Δ
t
,
S
h
)
+
k
(
∇
(
T
h
n
+
1
+
T
h
n
2
)
,
∇
S
h
)
+
μ
h
(
∇
T
h
n
+
1
,
∇
S
h
)
+
c
-
(
E
[
u
h
n
,
u
h
n
-
1
]
,
T
h
n
+
1
+
T
h
n
2
,
S
h
)
=
(
γ
(
t
n
+
1
/
2
)
,
S
h
)
+
μ
h
(
∇
T
h
n
,
∇
S
h
)
.We find that CNLES FE method requires the solution of only one linear problem at each time step; thus it needs less time than fully implicit Crank-Nicolson scheme. Denotek
min
⁡
=
min
⁡
(
k
e
,
k
s
), and k
max
⁡
=
max
⁡
(
k
e
,
k
s
).


## 3.2. Stability of the Method

Theorem 5.
Suppose0
<
k
min
⁡
≤
k
max
⁡
<
∞, and γ
∈
L
2
(
0
,
T
;
H
-
1
(
Ω
)
). Algorithm 4 is unconditionally stable in the following sense, for any h, Δ
t
>
0, μ
⩾
0, and 0
⩽
l
⩽
N
-
1:(20)
∥
T
h
l
+
1
∥
2
+
k
min
⁡
Δ
t
∑
n
=
0
l
∥
∇
(
T
h
n
+
1
+
T
h
n
2
)
∥
2
+
μ
h
Δ
t
∥
∇
T
h
l
+
1
∥
2
⩽
∥
T
h
0
∥
2
+
μ
h
Δ
t
∥
∇
T
h
0
∥
2
+
k
min
⁡
-
1
Δ
t
∑
n
=
0
l
∥
γ
(
t
n
+
1
/
2
)
∥
-
1
2
,
(21)
∥
u
h
l
+
1
∥
2
+
Pr
Δ
t
∑
n
=
0
l
∥
∇
(
u
h
n
+
1
+
u
h
n
2
)
∥
2
+
μ
h
Δ
t
∥
∇
u
h
l
+
1
∥
2
⩽
∥
u
h
0
∥
2
+
μ
h
Δ
t
∥
∇
u
h
0
∥
2
+
Pr
Ra
2
k
min
⁡
-
1
C
P
[
∑
n
=
0
l
k
min
⁡
-
1
∥
γ
(
t
n
+
(
1
/
2
)
)
∥
-
1
2
∥
T
h
0
∥
2
+
μ
h
Δ
t
∥
∇
T
h
0
∥
2
0000000000000000000
+
Δ
t
∑
n
=
0
l
k
min
⁡
-
1
∥
γ
(
t
n
+
1
/
2
)
∥
-
1
2
]
.Proof.
We first derive the stability of temperatureT and then of velocity u. Choosing S
h
=
(
T
h
1
+
T
h
0
)
/
2
∈
W
h in (17) yields
(22)
(
T
h
1
-
T
h
0
Δ
t
,
T
h
1
+
T
h
0
2
)
+
k
∥
∇
(
T
h
1
+
T
h
0
2
)
∥
2
+
μ
h
Δ
t
(
∇
(
T
h
1
-
T
h
0
Δ
t
)
,
∇
(
T
h
1
+
T
h
0
2
)
)
=
(
γ
(
t
1
/
2
)
,
T
h
1
+
T
h
0
2
)
.
Applying the Cauchy-Schwarz and Young’s inequalities gives
(23)
∥
T
h
1
∥
2
-
∥
T
h
0
∥
2
2
Δ
t
+
k
min
⁡
∥
∇
(
T
h
1
+
T
h
0
2
)
∥
2
+
μ
h
Δ
t
∥
∇
T
h
1
∥
2
-
∥
T
h
0
∥
2
2
Δ
t
⩽
1
2
k
min
⁡
-
1
∥
γ
(
t
1
/
2
)
∥
-
1
2
+
1
2
k
min
⁡
∥
∇
(
T
h
1
+
T
h
0
2
)
∥
2
.
Thus, we get the bound of temperature on the first time level:
(24)
∥
T
h
1
∥
2
+
Δ
t
k
min
⁡
∥
∇
(
T
h
1
+
T
h
0
2
)
∥
2
+
μ
h
Δ
t
∥
∇
T
h
1
∥
2
⩽
∥
T
h
0
∥
2
+
μ
h
Δ
t
∥
∇
T
h
0
∥
2
+
k
min
⁡
-
1
Δ
t
∥
γ
(
t
1
/
2
)
∥
-
1
2
.
For n
⩾
1, taking S
h
=
(
T
h
n
+
1
+
T
h
n
)
/
2
∈
W
h in (19) and using the Cauchy-Schwarz and Young inequalities lead to
(25)
(
∥
T
h
n
+
1
∥
2
-
∥
T
h
n
∥
2
)
+
k
min
⁡
Δ
t
∥
∇
(
T
h
n
+
1
+
T
h
n
2
)
∥
2
+
μ
h
Δ
t
∥
∇
T
h
n
+
1
∥
2
⩽
μ
h
Δ
t
∥
∇
T
h
n
∥
2
+
k
min
⁡
-
1
Δ
t
∥
γ
(
t
n
+
1
/
2
)
∥
-
1
2
.
Summing the above equation over n from 1 to l yields
(26)
∥
T
h
l
+
1
∥
2
+
k
min
⁡
Δ
t
∑
n
=
1
l
∥
∇
(
T
h
n
+
1
+
T
h
n
2
)
∥
2
+
μ
h
Δ
t
∥
∇
T
h
l
+
1
∥
2
⩽
∥
T
h
1
∥
2
+
μ
h
Δ
t
∥
∇
T
h
1
∥
2
+
k
min
⁡
-
1
Δ
t
∑
n
=
1
l
∥
γ
(
t
n
+
1
/
2
)
∥
-
1
2
.
Substituting the bound ∥
T
h
1
∥
2
+
Δ
t
k
min
⁡
∥
∇
(
(
T
h
1
+
T
h
0
)
/
2
)
∥
2
+
μ
h
Δ
t
∥
∇
T
h
1
∥
2 of (24) into the above relation, we obtain the result of (20).
Now we derive the stability of velocityu. Taking v
h
=
(
u
h
1
+
u
h
0
)
/
2
∈
V
h in (16) and applying the Cauchy-Schwarz and Young inequalities, we have
(27)
∥
u
h
1
∥
2
-
∥
u
h
0
∥
2
2
Δ
t
+
Pr
∥
∇
(
u
h
1
+
u
h
0
2
)
∥
2
+
μ
h
Δ
t
∥
∇
u
h
1
∥
2
-
∥
∇
u
h
0
∥
2
2
Δ
t
⩽
Pr
2
∥
∇
(
u
h
1
+
u
h
0
2
)
∥
2
+
Pr
R
a
2
2
∥
T
h
1
+
T
h
0
2
∥
-
1
2
.
Then
(28)
∥
u
h
1
∥
2
+
Pr
Δ
t
∥
∇
(
u
h
1
+
u
h
0
2
)
∥
2
+
μ
h
Δ
t
∥
∇
u
h
1
∥
2
⩽
∥
u
h
0
∥
2
+
μ
h
Δ
t
∥
∇
u
h
0
∥
2
+
Pr
R
a
2
Δ
t
∥
T
h
1
+
T
h
0
2
∥
-
1
2
.
For n
⩾
1, taking v
h
=
(
u
h
n
+
1
+
u
h
n
)
/
2
∈
V
h in (18) and making use of the Cauchy-Schwarz inequality and the Young’s inequality lead to
(29)
∥
u
h
n
+
1
∥
2
-
∥
u
h
n
∥
2
2
Δ
t
+
Pr
∥
∇
(
u
h
n
+
1
+
u
h
n
2
)
∥
2
+
μ
h
Δ
t
∥
∇
u
h
n
+
1
∥
2
-
∥
∇
u
h
n
∥
2
2
Δ
t
⩽
Pr
2
∥
∇
(
u
h
n
+
1
+
u
h
n
2
)
∥
2
+
Pr
R
a
2
2
∥
T
h
n
+
1
+
T
h
n
2
∥
-
1
2
.
Multiplying the above relation by 2
Δ
t and summing over n from 1 to l give
(30)
∥
u
h
l
+
1
∥
2
+
Pr
Δ
t
∑
n
=
1
l
∥
∇
(
u
h
n
+
1
+
u
h
n
2
)
∥
2
+
μ
h
Δ
t
∥
∇
u
h
l
+
1
∥
2
⩽
∥
u
h
1
∥
2
+
μ
h
Δ
t
∥
∇
u
h
1
∥
2
+
Δ
t
Pr
Ra
2
∑
n
=
1
l
∥
T
h
n
+
1
+
T
h
n
2
∥
-
1
2
.
Inserting the bound ∥
u
h
1
∥
2
+
Pr
Δ
t
∥
∇
(
(
u
h
1
+
u
h
0
)
/
2
)
∥
2
+
μ
h
Δ
t
∥
∇
u
h
1
∥
2 of (28) into the above relation and using the result (20), we get the desired result (21).


## 4. Error Estimate

Denoting(
u
(
t
n
+
1
)
,
p
(
t
n
+
1
)
,
T
(
t
n
+
1
)
)
=
(
u
n
+
1
,
p
n
+
1
,
T
n
+
1
) and splitting the error terms e
u
n
+
1, e
p
n
+
1, and e
T
n
+
1 into the approximation errors η
u
n
+
1, η
p
n
+
1, and η
T
n
+
1 and the model errors φ
u
n
+
1, φ
p
n
+
1, and φ
T
n
+
1, respectively,
(31)
e
u
n
+
1
=
u
n
+
1
-
u
h
n
+
1
=
(
u
n
+
1
-
u
~
n
+
1
)
-
(
u
h
n
+
1
-
u
~
n
+
1
)
=
η
u
n
+
1
-
φ
u
n
+
1
,
e
p
n
+
1
=
p
n
+
1
-
p
h
n
+
1
=
(
p
n
+
1
-
p
~
n
+
1
)
-
(
p
h
n
+
1
-
p
~
n
+
1
)
=
η
p
n
+
1
-
φ
p
n
+
1
,
e
T
n
+
1
=
T
n
+
1
-
T
h
n
+
1
=
(
T
n
+
1
-
T
~
n
+
1
)
-
(
T
h
n
+
1
-
T
~
n
+
1
)
=
η
T
n
+
1
-
φ
T
n
+
1
,
where (
u
~
n
+
1
,
p
~
n
+
1
,
T
~
n
+
1
) are the modified Stokes projections of (
u
(
t
n
+
1
)
,
p
(
t
n
+
1
)
,
T
(
t
n
+
1
)
) into spaces (X
h
,
M
h
,
W
h), respectively. For m
⩾
2, we assume the exact solution of problem (1) satisfies following regularity assumptions:
(32)
u
∈
L
2
(
0
,
t
*
;
H
m
+
1
(
Ω
)
d
)
∩
L
∞
(
0
,
t
*
;
H
m
+
1
(
Ω
)
d
)
,
u
t
∈
L
2
(
0
,
t
*
;
H
m
+
1
(
Ω
)
d
)
,
u
t
t
∈
L
2
(
0
,
t
*
;
H
2
(
Ω
)
d
)
,
u
t
t
t
∈
L
∞
(
0
,
t
*
;
L
2
(
Ω
)
d
)
,
p
t
t
∈
L
2
(
0
,
t
*
;
L
2
(
Ω
)
)
,
T
∈
L
2
(
0
,
t
*
;
H
m
+
1
(
Ω
)
)
∩
L
∞
(
0
,
t
*
;
H
m
+
1
(
Ω
)
)
,
T
t
∈
L
2
(
0
,
t
*
;
H
m
+
1
(
Ω
)
)
,
T
t
t
∈
L
2
(
0
,
t
*
;
H
2
(
Ω
)
)
,
T
t
t
t
∈
L
∞
(
0
,
t
*
;
L
2
(
Ω
)
)
.Theorem 6.
Assume that the solution(
u
,
p
,
T
) of problem (1) satisfies regularity assumptions (32). Let the finite element spaces (
X
h
,
M
h
,
W
h
) include continuous piecewise polynomials of degree m, m
-
1, and m (m
⩾
2), respectively. If Δ
t satisfies the condition
(33)
max
⁡
⁡
{
C
h
m
-
d
/
2
Δ
t
∥
u
∥
L
∞
(
0
,
t
*
;
H
m
+
1
(
Ω
)
)
,
00000
6
P
r
R
a
2
C
P
2
Δ
t
+
C
h
m
-
d
/
2
Δ
t
∥
T
∥
L
∞
(
0
,
t
*
;
H
m
+
1
(
Ω
)
)
}
⩽
1
2
,
then there is a constant C
^
=
C
^
(
P
r
,
k
,
R
a
,
Ω
,
m
,
u
,
p
,
T
) such that, for any 0
⩽
l
⩽
N
-
1,
(34)
∥
u
(
t
l
+
1
)
-
u
h
l
+
1
∥
+
(
Δ
t
∑
n
=
0
l
P
r
∥
∇
(
(
(
u
(
t
n
+
1
)
-
u
h
n
+
1
)
+
(
u
(
t
n
)
-
u
h
n
)
)
v
v
v
v
v
v
v
v
v
v
v
v
w
v
v
v
v
v
v
×
(
2
)
-
1
(
(
u
(
t
n
+
1
)
-
u
h
n
+
1
)
+
(
u
(
t
n
)
-
u
h
n
)
)
)
∥
2
∑
n
=
0
l
Pr
∥
∇
(
(
(
u
(
t
n
+
1
)
-
u
h
n
+
1
)
+
(
u
(
t
n
)
-
u
h
n
)
)
)
1
/
2
+
∥
T
(
t
l
+
1
)
-
T
h
l
+
1
∥
+
(
Δ
t
∑
n
=
0
l
k
min
⁡
∥
∇
(
(
(
T
(
t
n
+
1
)
-
T
h
n
+
1
)
+
(
T
(
t
n
)
-
T
h
n
)
)
v
v
v
v
v
v
v
v
v
v
v
v
w
v
v
v
v
v
v
v
×
(
2
)
-
1
(
(
T
(
t
n
+
1
)
-
T
h
n
+
1
)
+
(
T
(
t
n
)
-
T
h
n
)
)
)
∥
2
∑
n
=
0
l
)
1
/
2
+
(
μ
h
Δ
t
)
1
/
2
∥
∇
(
u
(
t
l
+
1
)
-
u
h
l
+
1
)
∥
+
(
μ
h
Δ
t
)
1
/
2
∥
∇
(
T
(
t
l
+
1
)
-
T
h
l
+
1
)
∥
⩽
C
^
(
h
m
+
μ
h
Δ
t
+
Δ
t
2
)
.Remark 7.
We see that our proof is conditionally convergent. We think that the condition (33) is not necessary theoretical, and it might be removed. But we cannot make it herein and will study it in the future.Now we give the outline of the proof. The proof is proved in Part 1 (establishing the error estimate for Step 1 in Algorithm4) and Part 2 (deriving the error estimate for Step 2 in Algorithm 4). In each parts, we first derive the error equations of the momentum equation and the energy equation; then give the error estimates of the momentum error equation and the energy error equation, respectively, and finally derive Theorem 6 by using the Gronwall lemma and the triangle inequality.Proof.
Part
1. Derive the error estimate of Step 1 in Algorithm 4.
Step
1.1. Establishe the error equations of the momentum equation and the energy equation. Subtracting equations (16) and (17) from (4) at t
=
t
1
/
2 for ∀
(
v
h
,
S
h
)
∈
(
V
h
,
W
h
) gives
(35)
(
u
t
(
t
1
/
2
)
-
u
h
1
-
u
h
0
Δ
t
,
v
h
)
+
Pr
(
∇
u
(
t
1
/
2
)
-
∇
(
u
h
1
+
u
h
0
2
)
,
∇
v
h
)
-
μ
h
(
∇
(
u
h
1
-
u
h
0
)
,
∇
v
h
)
+
c
(
u
(
t
1
/
2
)
,
u
(
t
1
/
2
)
,
v
h
)
-
c
(
u
h
0
,
u
h
1
+
u
h
0
2
,
v
h
)
-
(
p
(
t
1
/
2
)
-
p
h
1
+
p
h
n
2
,
∇
·
v
h
)
=
Pr
Ra
(
ζ
T
(
t
1
/
2
)
-
ζ
T
h
1
+
T
h
0
2
,
v
h
)
,
(36)
(
T
t
(
t
1
/
2
)
-
T
h
1
-
T
h
0
Δ
t
,
S
h
)
+
k
(
∇
T
(
t
1
/
2
)
-
∇
(
T
h
1
+
T
h
0
2
)
,
∇
S
h
)
-
μ
h
(
∇
(
T
h
1
-
T
h
0
)
,
∇
S
h
)
+
c
-
(
u
(
t
1
/
2
)
,
T
(
t
1
/
2
)
,
S
h
)
-
c
-
(
u
h
0
,
T
h
1
+
T
h
0
2
,
S
h
)
=
0
.
Adding and subtracting
(37)
(
u
(
t
1
)
-
u
(
t
0
)
Δ
t
,
v
h
)
+
Pr
(
∇
(
u
(
t
1
)
+
u
(
t
0
)
2
)
,
∇
v
h
)
+
μ
h
(
∇
u
(
t
1
)
-
∇
u
(
t
0
)
,
∇
v
h
)
+
(
p
(
t
1
)
+
p
(
t
0
)
2
,
∇
·
v
h
)
-
Pr
Ra
(
ζ
T
(
t
1
)
+
T
(
t
0
)
2
,
v
h
)
,
(
T
(
t
1
)
-
T
(
t
0
)
Δ
t
,
S
h
)
+
k
(
∇
(
T
(
t
1
)
+
T
(
t
0
)
2
)
,
∇
S
h
)
+
μ
h
(
∇
(
T
(
t
1
)
-
T
(
t
0
)
)
,
∇
S
h
)
to (35) and (36), respectively, we have
(38)
(
e
u
1
-
e
u
0
Δ
t
,
v
h
)
+
Pr
(
∇
e
u
1
/
2
,
∇
v
h
)
+
μ
h
(
∇
e
u
1
-
∇
e
u
0
,
∇
v
h
)
=
(
e
p
1
/
2
,
∇
·
v
h
)
+
M
O
1
(
u
,
p
,
T
,
e
T
;
v
h
)
,
(
e
T
1
-
e
T
0
Δ
t
,
S
h
)
+
k
(
∇
e
T
1
/
2
,
∇
S
h
)
+
μ
h
(
∇
(
e
T
1
-
e
T
0
)
,
∇
S
h
)
=
E
N
1
(
u
,
T
;
S
h
)
,
where
(39)
M
O
1
(
u
,
p
,
T
,
e
T
;
v
h
)
=
(
u
(
t
1
)
-
u
(
t
0
)
Δ
t
-
u
t
(
t
1
/
2
)
,
v
h
)
+
Pr
(
∇
(
u
(
t
1
)
+
u
(
t
0
)
2
-
u
(
t
1
/
2
)
)
,
∇
v
h
)
+
μ
h
(
∇
u
(
t
1
)
-
∇
u
(
t
0
)
,
∇
v
h
)
+
(
p
(
t
1
/
2
)
-
p
(
t
1
)
+
p
(
t
0
)
2
,
∇
·
v
h
)
-
Pr
Ra
(
ζ
T
(
t
1
)
+
T
(
t
0
)
2
,
v
h
)
+
Pr
Ra
(
ζ
T
(
t
1
/
2
)
,
v
h
)
+
Pr
Ra
(
ζ
e
T
1
/
2
,
v
h
)
-
c
(
u
(
t
1
/
2
)
,
u
(
t
1
/
2
)
,
v
h
)
+
c
(
u
h
0
,
u
h
1
+
u
h
0
2
,
v
h
)
,
E
N
1
(
u
,
T
;
S
h
)
=
(
T
(
t
1
)
-
T
(
t
0
)
Δ
t
-
T
t
(
t
1
/
2
)
,
S
h
)
+
k
(
∇
(
T
(
t
1
)
+
T
(
t
0
)
2
)
-
∇
T
(
t
1
/
2
)
,
∇
S
h
)
+
μ
h
(
∇
T
(
t
1
)
-
∇
T
(
t
0
)
,
∇
S
h
)
-
c
-
(
u
(
t
1
/
2
)
,
T
(
t
1
/
2
)
,
S
h
)
+
c
-
(
u
h
0
,
T
h
1
+
T
h
0
2
,
S
h
)
.
Using the error decomposition (31) and choosing v
h
=
φ
u
1
/
2, S
h
=
φ
T
1
/
2 give
(40)
1
2
Δ
t
(
∥
φ
u
1
∥
2
-
∥
φ
u
0
∥
2
)
+
Pr
∥
∇
φ
u
1
/
2
∥
2
+
μ
h
2
(
∥
∇
φ
u
1
∥
2
-
∥
∇
φ
u
0
∥
2
)
=
(
η
u
1
-
η
u
0
Δ
t
,
φ
u
1
/
2
)
+
Pr
(
∇
η
u
1
/
2
,
∇
φ
u
1
/
2
)
+
μ
h
Δ
t
(
∇
(
η
u
1
-
η
u
0
Δ
t
)
,
∇
φ
u
1
/
2
)
-
(
e
p
1
/
2
,
∇
·
φ
u
1
/
2
)
-
M
O
1
(
u
,
p
,
T
,
e
T
;
φ
u
1
/
2
)
,
(41)
1
2
Δ
t
(
∥
φ
T
1
∥
2
-
∥
φ
T
0
∥
2
)
+
k
∥
∇
φ
T
1
/
2
∥
2
+
μ
h
2
(
∥
∇
φ
T
1
∥
2
-
∥
∇
φ
T
0
∥
2
)
=
(
η
T
1
-
η
T
0
Δ
t
,
φ
T
1
/
2
)
+
k
(
∇
η
T
1
/
2
,
∇
φ
T
1
/
2
)
+
μ
h
Δ
t
(
∇
(
η
T
1
-
η
T
0
Δ
t
)
,
∇
φ
T
1
/
2
)
-
E
N
1
(
u
,
T
;
φ
T
1
/
2
)
.
Step
1.2. Derive the error inequalities of the momentum equation (40). From the definition of the modified Stokes projection (11)–(13), we have
(42)
Pr
(
∇
η
u
1
/
2
,
∇
φ
u
1
/
2
)
-
(
e
p
1
/
2
,
∇
·
φ
u
1
/
2
)
=
0
,
k
(
∇
η
T
1
/
2
,
∇
φ
T
1
/
2
)
=
0
.
For the linear terms on the RHS of (40), applying the Cauchy-Schwarz, Young inequalities. and the approximation properties (9) gives
(43)
|
(
η
u
1
-
η
u
0
Δ
t
,
φ
u
1
/
2
)
|
⩽
Pr
24
∥
∇
φ
u
1
/
2
∥
2
+
6
C
P
2
Pr
∥
η
u
1
-
η
u
0
Δ
t
∥
2
⩽
Pr
24
∥
∇
φ
u
1
/
2
∥
2
+
C
h
2
m
+
2
Pr
Δ
t
∥
u
t
∥
L
2
(
t
0
,
t
1
;
H
m
+
1
(
Ω
)
)
2
,
|
μ
h
Δ
t
(
∇
(
η
u
1
-
η
u
0
Δ
t
)
,
∇
φ
u
1
/
2
)
|
⩽
Pr
24
∥
∇
φ
u
1
/
2
∥
2
+
6
(
μ
h
Δ
t
)
2
Pr
∥
∇
(
η
u
1
-
η
u
0
Δ
t
)
∥
2
⩽
Pr
24
∥
∇
φ
u
1
/
2
∥
2
+
C
μ
2
h
2
m
+
2
Δ
t
Pr
∥
u
t
∥
L
2
(
t
0
,
t
1
;
H
m
+
1
(
Ω
)
)
2
.
We next estimate the term M
O
1
(
u
,
p
,
T
,
e
T
;
φ
u
1
/
2
) of (40). For the linear terms of M
O
1
(
u
,
p
,
T
,
e
T
;
φ
u
1
/
2
), using the Cauchy-Schwarz and Young’s inequalities together with Taylor expansion on t, we get
(44)
|
(
u
(
t
1
)
-
u
(
t
0
)
Δ
t
-
u
t
(
t
1
/
2
)
,
φ
u
1
/
2
)
|
⩽
Pr
24
∥
∇
φ
u
1
/
2
∥
2
+
C
Δ
t
4
Pr
max
⁡
t
∈
[
t
0
,
t
1
]
∥
u
t
t
t
(
t
)
∥
2
,
Pr
|
(
∇
(
u
(
t
1
)
+
u
(
t
0
)
2
-
u
(
t
1
/
2
)
)
,
∇
φ
u
1
/
2
)
|
⩽
Pr
24
∥
∇
φ
u
1
/
2
∥
2
+
C
Pr
Δ
t
4
max
⁡
t
∈
[
t
0
,
t
1
]
∥
∇
u
t
t
(
t
)
∥
2
,
μ
h
Δ
t
|
(
∇
(
u
(
t
1
)
-
u
(
t
0
)
Δ
t
)
,
∇
φ
u
1
/
2
)
|
⩽
Pr
24
∥
∇
φ
u
1
/
2
∥
2
+
C
μ
2
h
2
Δ
t
2
Pr
max
⁡
t
∈
[
t
0
,
t
1
]
∥
∇
u
t
(
t
)
∥
2
,
|
(
p
(
t
1
)
+
p
(
t
0
)
2
-
p
(
t
1
/
2
)
,
∇
·
φ
u
1
/
2
)
|
⩽
Pr
24
∥
∇
φ
u
1
/
2
∥
2
+
C
Δ
t
4
Pr
max
⁡
t
∈
[
t
0
,
t
1
]
∥
p
t
t
(
t
)
∥
2
,
Pr
Ra
|
(
ζ
T
(
t
1
)
+
T
(
t
0
)
2
-
ζ
T
(
t
1
/
2
)
,
φ
u
1
/
2
)
|
⩽
Pr
24
∥
∇
φ
u
1
/
2
∥
2
+
C
Pr
Ra
2
Δ
t
4
max
⁡
t
∈
[
t
0
,
t
1
]
∥
T
t
t
(
t
)
∥
2
,
Pr
Ra
|
(
ζ
e
T
1
/
2
,
φ
u
1
/
2
)
|
=
Pr
Ra
|
(
ζ
η
T
1
/
2
-
ζ
φ
T
1
/
2
,
φ
u
1
/
2
)
|
⩽
Pr
24
∥
∇
φ
u
1
/
2
∥
2
+
6
Pr
Ra
2
C
P
2
∥
η
T
1
/
2
∥
2
+
Pr
24
∥
∇
φ
u
1
/
2
∥
2
+
6
4
Pr
Ra
2
C
P
2
∥
φ
T
1
∥
2
⩽
Pr
12
∥
∇
φ
u
1
/
2
∥
2
+
6
Pr
Ra
2
C
P
2
∥
η
T
1
/
2
∥
2
+
1
8
Δ
t
∥
φ
T
1
∥
2
,
where we have used the assumption (33) in the last equation. For the the nonlinear terms of |
M
O
1
(
u
,
p
,
T
,
e
T
;
φ
u
1
/
2
)
| in (40), by adding and subtracting some terms and using Lemma 1, we have
(45)
-
c
(
u
(
t
1
/
2
)
,
u
(
t
1
/
2
)
,
φ
u
1
/
2
)
+
c
(
u
h
0
,
u
h
1
+
u
h
0
2
,
φ
u
1
/
2
)
=
-
c
(
u
(
t
1
/
2
)
,
u
(
t
1
/
2
)
,
φ
u
1
/
2
)
-
c
(
u
h
0
,
e
u
1
/
2
,
φ
u
1
/
2
)
-
c
(
u
(
t
0
)
-
u
h
0
,
u
(
t
1
)
+
u
(
t
0
)
2
,
φ
u
1
/
2
)
+
c
(
u
(
t
0
)
,
u
(
t
1
)
+
u
(
t
0
)
2
,
φ
u
1
/
2
)
.
Using the inequalities of Lemma 1 and the regularity assumptions (32) on u, then the second and the third terms of (45) are bounded by
(46)
|
c
(
u
h
0
,
e
u
1
/
2
,
φ
u
1
/
2
)
|
=
|
c
(
u
h
0
,
η
u
1
/
2
,
φ
u
1
/
2
)
|
=
|
-
c
(
η
u
0
,
η
u
1
/
2
,
φ
u
1
/
2
)
+
c
(
u
(
t
0
)
,
η
u
1
/
2
,
φ
u
1
/
2
)
|
⩽
Pr
24
∥
∇
φ
u
1
/
2
∥
2
+
C
P
r
-
1
∥
∇
η
u
1
/
2
∥
2
∥
∇
η
u
0
∥
2
+
C
P
r
-
1
∥
∇
η
u
1
/
2
∥
2
,
|
c
(
u
(
t
0
)
-
u
h
0
,
u
(
t
1
)
+
u
(
t
0
)
2
,
φ
u
1
/
2
)
|
=
|
c
(
η
u
0
,
u
(
t
1
)
+
u
(
t
0
)
2
,
φ
u
1
/
2
)
|
⩽
C
∥
∇
φ
u
1
/
2
∥
∥
∇
η
u
0
∥
⩽
Pr
24
∥
∇
φ
u
1
/
2
∥
2
+
C
P
r
-
1
∥
∇
η
u
0
∥
2
.
Using Taylor expansion on t, we write the first and the fourth terms of (45) as
(47)
|
(
u
(
t
0
)
,
u
(
t
1
)
+
u
(
t
0
)
2
,
φ
u
1
/
2
)
-
c
(
u
(
t
1
/
2
)
,
u
(
t
1
/
2
)
,
φ
u
1
/
2
)
00
+
c
(
u
(
t
0
)
,
u
(
t
1
)
+
u
(
t
0
)
2
,
φ
u
1
/
2
)
|
⩽
|
c
(
u
(
t
0
)
-
u
(
t
1
/
2
)
,
u
(
t
1
/
2
)
,
φ
u
1
/
2
)
|
+
C
Δ
t
2
|
c
(
u
(
t
0
)
,
u
t
t
(
t
θ
)
,
φ
u
1
/
2
)
|
⩽
C
Δ
t
|
c
(
u
t
(
t
θ
)
,
u
(
t
1
/
2
)
,
φ
u
1
/
2
)
|
+
C
Δ
t
2
|
c
(
u
(
t
0
)
,
u
t
t
(
t
θ
)
,
φ
u
1
/
2
)
|
with t
θ
∈
(
t
0
,
t
1
). Applying Lemma 1, Taylor expansion on t, and the regularity assumptions (32), the terms of (47) are bounded by
(48)
Δ
t
|
c
(
u
t
(
t
θ
)
,
u
(
t
1
/
2
)
,
φ
u
1
/
2
)
|
⩽
C
Δ
t
(
∥
u
t
(
t
θ
)
∥
∥
∇
u
(
t
1
/
2
)
∥
L
∞
(
Ω
)
00000000
+
∥
∇
u
t
(
t
θ
)
∥
∥
u
(
t
1
/
2
)
∥
L
∞
(
Ω
)
)
∥
φ
u
1
/
2
∥
⩽
1
4
Δ
t
∥
φ
u
1
∥
2
+
C
Δ
t
3
(
∥
u
t
(
t
θ
)
∥
∥
∇
u
(
t
1
/
2
)
∥
L
∞
(
Ω
)
0000000000
+
∥
∇
u
t
(
t
θ
)
∥
∥
u
(
t
1
/
2
)
∥
L
∞
(
Ω
)
)
2
,
C
Δ
t
2
|
c
(
u
(
t
0
)
,
u
t
t
(
t
θ
)
,
φ
u
1
/
2
)
|
⩽
Pr
24
∥
∇
φ
u
1
/
2
∥
2
+
C
Δ
t
4
Pr
∥
∇
u
(
t
0
)
∥
2
∥
∇
u
t
t
(
t
θ
)
∥
2
.
Combining these inequalities with (40), we obtain
(49)
1
2
Δ
t
(
∥
φ
u
1
∥
2
-
∥
φ
u
0
∥
2
)
+
Pr
2
∥
∇
φ
u
1
/
2
∥
2
vvvvv
+
μ
h
2
(
∥
∇
φ
u
1
∥
2
-
∥
∇
φ
u
0
∥
2
)
⩽
1
4
Δ
t
∥
φ
u
1
∥
2
+
1
8
Δ
t
∥
φ
T
1
∥
2
+
C
h
2
m
+
2
Pr
(
1
Δ
t
+
μ
2
Δ
t
)
∥
u
t
∥
L
2
(
t
0
,
t
1
;
H
m
+
1
(
Ω
)
)
2
+
C
P
r
-
1
∥
∇
η
u
1
/
2
∥
2
∥
∇
η
u
0
∥
2
+
C
P
r
-
1
∥
∇
η
u
1
/
2
∥
2
+
C
P
r
-
1
∥
∇
η
u
0
∥
2
+
C
Δ
t
4
P
r
-
1
[
max
⁡
t
∈
[
t
0
,
t
1
]
∥
u
t
t
t
(
t
)
∥
2
0000000000000000
+
P
r
2
max
⁡
t
∈
[
t
0
,
t
1
]
∥
∇
u
t
t
(
t
)
∥
2
0000000000000000
+
max
⁡
t
∈
[
t
0
,
t
1
]
∥
p
t
t
(
t
)
∥
2
0000000000000000
+
P
r
2
R
a
2
max
⁡
t
∈
[
t
0
,
t
1
]
∥
T
t
t
(
t
)
∥
2
0000000000000000
+
∥
∇
u
(
t
0
)
∥
2
∥
∇
u
t
t
(
t
θ
)
∥
2
max
⁡
t
∈
[
t
0
,
t
1
]
]
+
6
Pr
R
a
2
C
P
2
∥
η
T
1
/
2
∥
2
+
C
P
r
-
1
μ
2
h
2
Δ
t
2
max
⁡
t
∈
[
t
0
,
t
1
]
∥
∇
u
t
(
t
)
∥
2
+
C
Δ
t
3
(
∥
u
t
(
t
θ
)
∥
∥
∇
u
(
t
1
/
2
)
∥
L
∞
(
Ω
)
vvvvvvvvvvv
+
∥
∇
u
t
(
t
θ
)
∥
∥
u
(
t
1
/
2
)
∥
L
∞
(
Ω
)
)
2
.
Step
1.3. Derive the error inequality of the energy error equation (41). Similarly, for the linear terms on the RHS of (41), we have
(50)
|
(
η
T
1
-
η
T
0
Δ
t
,
φ
T
1
/
2
)
|
⩽
k
min
⁡
16
∥
∇
φ
T
1
/
2
∥
2
+
4
k
min
⁡
-
1
C
P
∥
η
T
1
-
η
T
0
Δ
t
∥
2
⩽
k
min
⁡
16
∥
∇
φ
T
1
/
2
∥
2
+
C
h
2
m
+
2
k
min
⁡
Δ
t
∥
T
t
∥
L
2
(
t
0
,
t
1
;
H
m
+
1
(
Ω
)
)
2
,
|
μ
h
Δ
t
(
∇
(
η
T
1
-
η
T
0
Δ
t
)
,
∇
φ
T
1
/
2
)
|
⩽
k
min
⁡
16
∥
∇
φ
T
1
/
2
∥
2
+
4
(
μ
h
Δ
t
)
2
k
min
⁡
∥
∇
(
η
T
1
-
η
T
0
Δ
t
)
∥
2
⩽
k
min
⁡
16
∥
∇
φ
T
1
/
2
∥
2
+
C
μ
2
h
2
m
+
2
Δ
t
k
min
⁡
∥
T
t
∥
L
2
(
t
0
,
t
1
;
H
m
+
1
(
Ω
)
)
2
.
For the linear terms of E
N
1
(
u
,
T
;
φ
T
1
/
2
), we arrive at
(51)
|
(
T
(
t
1
)
-
T
(
t
0
)
Δ
t
-
T
t
(
t
1
/
2
)
,
φ
T
1
/
2
)
|
⩽
k
min
⁡
16
∥
∇
φ
T
1
/
2
∥
2
+
C
Δ
t
4
k
min
⁡
max
⁡
t
∈
[
t
0
,
t
1
]
∥
T
t
t
t
(
t
)
∥
2
,
k
|
(
∇
(
T
(
t
1
)
+
T
(
t
0
)
2
-
T
(
t
1
/
2
)
)
,
∇
φ
T
1
/
2
)
|
⩽
k
min
⁡
16
∥
∇
φ
T
1
/
2
∥
2
+
C
k
max
⁡
2
k
min
⁡
Δ
t
4
max
⁡
t
∈
[
t
0
,
t
1
]
∥
∇
T
t
t
(
t
)
∥
2
,
μ
h
Δ
t
|
(
∇
(
T
(
t
1
)
-
T
(
t
0
)
Δ
t
)
,
∇
φ
T
1
/
2
)
|
⩽
k
min
⁡
16
∥
∇
φ
T
1
/
2
∥
2
+
C
μ
2
h
2
Δ
t
2
k
min
⁡
max
⁡
t
∈
[
t
0
,
t
1
]
∥
∇
T
t
(
t
)
∥
2
.
Using Lemma 1, we write the nonlinear terms of |
E
N
1
(
u
,
T
,
φ
T
1
/
2
)
| in (41) as
(52)
-
c
-
(
u
(
t
1
/
2
)
,
T
(
t
1
/
2
)
,
φ
T
1
/
2
)
+
c
-
(
u
h
0
,
T
h
1
+
T
h
0
2
,
φ
T
1
/
2
)
=
-
c
-
(
u
(
t
1
/
2
)
,
T
(
t
1
/
2
)
,
φ
T
1
/
2
)
-
c
-
(
u
h
0
,
e
T
1
/
2
,
φ
T
1
/
2
)
-
c
-
(
u
(
t
0
)
-
u
h
0
,
T
(
t
1
)
+
T
(
t
0
)
2
,
φ
T
1
/
2
)
+
c
-
(
u
(
t
0
)
,
T
(
t
1
)
+
T
(
t
0
)
2
,
φ
T
1
/
2
)
.
From Lemma 1, Taylor expansion on t, and the regularity assumptions (32), it follows that
(53)
|
c
-
(
u
h
0
,
e
T
1
/
2
,
φ
T
1
/
2
)
|
=
|
c
-
(
u
h
0
,
η
T
1
/
2
,
φ
T
1
/
2
)
|
=
|
-
c
-
(
η
u
0
,
η
T
1
/
2
,
φ
T
1
/
2
)
+
c
-
(
u
(
t
0
)
,
η
T
1
/
2
,
φ
T
1
/
2
)
|
⩽
k
min
⁡
16
∥
∇
φ
T
1
/
2
∥
2
+
C
k
min
⁡
-
1
∥
∇
η
T
1
/
2
∥
2
∥
∇
η
u
0
∥
2
+
C
k
min
⁡
-
1
∥
∇
η
T
1
/
2
∥
2
,
|
c
-
(
u
(
t
0
)
-
u
h
0
,
T
(
t
1
)
+
T
(
t
0
)
2
,
φ
T
1
/
2
)
|
=
|
c
-
(
η
u
0
,
T
(
t
1
)
+
T
(
t
0
)
2
,
φ
T
1
/
2
)
|
⩽
C
∥
∇
φ
T
1
/
2
∥
∥
∇
η
u
0
∥
⩽
k
min
⁡
16
∥
∇
φ
T
1
/
2
∥
2
+
C
k
min
⁡
-
1
∥
∇
η
u
0
∥
2
,
|
(
u
(
t
0
)
,
T
(
t
1
)
+
T
(
t
0
)
2
,
φ
T
1
/
2
)
-
c
-
(
u
(
t
1
/
2
)
,
T
(
t
1
/
2
)
,
φ
T
1
/
2
)
0
+
c
-
(
u
(
t
0
)
,
T
(
t
1
)
+
T
(
t
0
)
2
,
φ
T
1
/
2
)
|
⩽
1
8
Δ
t
∥
φ
T
1
∥
2
+
C
Δ
t
3
(
∥
u
t
(
t
θ
)
∥
∥
∇
T
(
t
1
/
2
)
∥
L
∞
(
Ω
)
vvvvvvvvvvvvv
+
∥
∇
u
t
(
t
θ
)
∥
∥
T
(
t
1
/
2
)
∥
L
∞
(
Ω
)
)
2
+
k
min
⁡
16
∥
∇
φ
T
1
/
2
∥
2
+
C
Δ
t
4
k
min
⁡
∥
∇
T
t
t
(
t
θ
)
∥
2
,
for any t
θ
∈
(
t
0
,
t
1
). Combining equations (52)-(53) with (41), we get
(54)
1
2
Δ
t
(
∥
φ
T
1
∥
2
-
∥
φ
T
0
∥
2
)
+
k
min
⁡
2
∥
∇
φ
T
1
/
2
∥
2
+
μ
h
2
(
∥
∇
φ
T
1
∥
2
-
∥
∇
φ
T
0
∥
2
)
⩽
1
8
Δ
t
∥
φ
T
1
∥
2
+
C
h
2
m
+
2
k
min
⁡
(
1
Δ
t
+
μ
2
Δ
t
)
∥
T
t
∥
L
2
(
t
0
,
t
1
;
H
m
+
1
(
Ω
)
)
2
+
C
k
min
⁡
(
max
⁡
t
∈
[
t
0
,
t
1
]
∥
∇
η
u
0
∥
2
+
∥
∇
η
T
1
/
2
∥
2
∥
∇
η
u
0
∥
2
000000000000
+
∥
∇
η
T
1
/
2
∥
2
+
μ
2
h
2
Δ
t
2
max
⁡
t
∈
[
t
0
,
t
1
]
∥
∇
T
t
(
t
)
∥
2
)
+
C
Δ
t
4
k
min
⁡
[
max
⁡
t
∈
[
t
0
,
t
1
]
∥
T
t
t
t
(
t
)
∥
2
+
k
max
⁡
2
max
⁡
t
∈
[
t
0
,
t
1
]
∥
∇
T
t
t
(
t
)
∥
2
0000000000000
+
∥
∇
T
t
t
(
t
θ
)
∥
2
max
⁡
t
∈
[
t
0
,
t
1
]
]
+
C
Δ
t
3
(
∥
u
t
(
t
θ
)
∥
∥
∇
T
(
t
1
/
2
)
∥
L
∞
(
Ω
)
000000000000
+
∥
∇
u
t
(
t
θ
)
∥
∥
T
(
t
1
/
2
)
∥
L
∞
(
Ω
)
)
2
.
Step
1.4. Combining (49) with (54), multiplying it by 2
Δ
t, and using the approximation properties (9) yield
(55)
1
2
∥
φ
u
1
∥
2
+
Pr
Δ
t
∥
∇
φ
u
1
/
2
∥
2
+
μ
h
Δ
t
(
∥
∇
φ
u
1
∥
2
-
∥
∇
φ
u
0
∥
2
)
+
1
2
∥
φ
T
1
∥
2
+
k
min
⁡
Δ
t
∥
∇
φ
T
1
/
2
∥
2
+
μ
h
Δ
t
(
∥
∇
φ
T
1
∥
2
-
∥
∇
φ
T
0
∥
2
)
⩽
C
h
2
m
+
2
(
1
+
μ
2
Δ
t
2
)
×
(
Pr
-
1
∥
u
t
∥
L
2
(
t
0
,
t
1
;
H
m
+
1
(
Ω
)
)
2
000000
+
k
min
⁡
-
1
∥
T
t
∥
L
2
(
t
0
,
t
1
;
H
m
+
1
(
Ω
)
)
2
)
000000
+
C
Δ
t
h
2
m
P
r
-
1
(
∥
u
(
t
1
)
∥
m
+
1
2
+
∥
u
(
t
0
)
∥
m
+
1
2
)
000000
+
C
Δ
t
h
2
m
k
min
⁡
-
1
(
∥
T
(
t
1
)
∥
m
+
1
2
+
∥
T
(
t
0
)
∥
m
+
1
2
0000000000000000000
+
∥
u
(
t
0
)
∥
m
+
1
2
)
+
C
Δ
t
5
P
r
-
1
[
max
⁡
t
∈
[
t
0
,
t
1
]
∥
u
t
t
t
(
t
)
∥
2
0000000000000000
+
P
r
2
max
⁡
t
∈
[
t
0
,
t
1
]
∥
∇
u
t
t
(
t
)
∥
2
0000000000000000
+
max
⁡
t
∈
[
t
0
,
t
1
]
∥
∇
p
t
t
(
t
)
∥
2
0000000000000000
+
P
r
2
R
a
2
max
⁡
t
∈
[
t
0
,
t
1
]
∥
T
t
t
(
t
)
∥
2
0000000000000000
+
∥
∇
u
(
t
0
)
∥
2
∥
∇
u
t
t
(
t
θ
)
∥
2
max
⁡
t
∈
[
t
0
,
t
1
]
]
+
C
Δ
t
5
k
min
⁡
-
1
[
max
⁡
t
∈
[
t
0
,
t
1
]
∥
T
t
t
t
(
t
)
∥
2
0000000000000000
+
k
max
⁡
2
max
⁡
t
∈
[
t
0
,
t
1
]
∥
∇
T
t
t
(
t
)
∥
2
0000000000000000
+
∥
∇
T
t
t
(
t
θ
)
∥
2
max
⁡
t
∈
[
t
0
,
t
1
]
]
+
C
Pr
R
a
2
C
P
2
Δ
t
h
2
m
+
2
(
∥
T
(
t
1
)
∥
m
+
1
2
+
∥
T
(
t
0
)
∥
m
+
1
2
)
+
C
μ
2
h
2
Δ
t
3
(
P
r
-
1
max
⁡
t
∈
[
t
0
,
t
1
]
∥
∇
u
t
(
t
)
∥
2
000000000000000
+
k
min
⁡
-
1
max
⁡
t
∈
[
t
0
,
t
1
]
∥
∇
T
t
(
t
)
∥
2
)
+
C
Δ
t
4
(
∥
u
t
(
t
θ
)
∥
∥
∇
u
(
t
1
/
2
)
∥
L
∞
(
Ω
)
00000000000
+
∥
∇
u
t
(
t
θ
)
∥
∥
u
(
t
1
/
2
)
∥
L
∞
(
Ω
)
)
2
+
C
Δ
t
4
(
∥
u
t
(
t
θ
)
∥
∥
∇
T
(
t
1
/
2
)
∥
L
∞
(
Ω
)
00000000000
+
∥
∇
u
t
(
t
θ
)
∥
∥
T
(
t
1
/
2
)
∥
L
∞
(
Ω
)
)
2
,
where we have used φ
u
1
/
2
=
(
1
/
2
)
φ
u
1, φ
T
1
/
2
=
(
1
/
2
)
φ
T
1. This ended the proof of Part 1.
Part
2. We derive the error estimate of Step 2 in Algorithm 4.
Step
2.1. Establis the error relations of Step 2. For n
⩾
1 and for all (
v
h
,
S
h
)
∈
(
V
h
,
W
h
), taking the variational formulation (1) at t
=
t
n
+
1
/
2 gives
(56)
(
u
t
(
t
n
+
1
/
2
)
,
v
h
)
+
Pr
(
∇
u
(
t
n
+
1
/
2
)
,
∇
v
h
)
+
c
(
u
(
t
n
+
1
/
2
)
,
u
(
t
n
+
1
/
2
)
,
v
h
)
-
(
p
(
t
n
+
1
/
2
)
,
∇
·
v
h
)
=
Pr
Ra
(
ζ
T
(
t
n
+
1
/
2
)
,
v
h
)
,
(
T
t
(
t
n
+
1
/
2
)
,
S
h
)
+
k
(
∇
T
(
t
n
+
1
/
2
)
,
∇
S
h
)
+
c
-
(
u
(
t
n
+
1
/
2
)
,
T
(
t
n
+
1
/
2
)
,
S
h
)
=
(
γ
(
t
n
+
1
/
2
)
,
S
h
)
.
Then, subtract (18)-(19) from (56) to get
(57)
(
u
t
(
t
n
+
1
/
2
)
-
u
h
n
+
1
-
u
h
n
Δ
t
,
v
h
)
+
Pr
⁡
(
∇
(
u
(
t
n
+
1
/
2
)
(
-
u
h
n
+
1
+
u
h
n
2
)
,
∇
v
h
)
)
-
μ
h
(
∇
(
u
h
n
+
1
-
u
h
n
)
,
∇
v
h
)
+
c
(
u
(
t
n
+
1
/
2
)
,
u
(
t
n
+
1
/
2
)
,
v
h
)
-
c
(
E
[
u
h
n
,
u
h
n
-
1
]
,
u
h
n
+
1
+
u
h
n
2
,
v
h
)
-
(
p
(
t
n
+
1
/
2
)
-
p
h
n
+
1
+
p
h
n
2
,
∇
·
v
h
)
=
Pr
Ra
(
ζ
T
(
t
n
+
1
/
2
)
-
ζ
T
h
n
+
1
+
T
h
n
2
,
v
h
)
,
(
T
t
(
t
n
+
1
/
2
)
-
T
h
n
+
1
-
T
h
n
Δ
t
,
S
h
)
+
k
(
∇
(
T
(
t
n
+
1
/
2
)
-
T
h
n
+
1
+
T
h
n
2
,
∇
S
h
)
)
-
μ
h
(
∇
(
T
h
n
+
1
-
T
h
n
)
,
∇
S
h
)
+
c
-
(
u
(
t
n
+
1
/
2
)
,
T
(
t
n
+
1
/
2
)
,
S
h
)
-
c
-
(
E
[
u
h
n
,
u
h
n
-
1
]
,
T
h
n
+
1
+
T
h
n
2
,
S
h
)
=
0
.
By adding and subtracting some terms to (57), we can obtain following error equations (recall that (
q
h
,
∇
·
v
h
)
=
0, ∀
q
h
∈
M
h):
(58)
1
Δ
t
(
e
u
n
+
1
-
e
u
n
,
v
h
)
+
Pr
(
∇
e
u
n
+
1
/
2
,
∇
v
h
)
+
μ
h
(
∇
e
u
n
+
1
-
∇
e
u
n
,
∇
v
h
)
=
(
e
p
n
+
1
/
2
,
∇
·
v
h
)
+
M
O
n
(
u
,
p
,
T
,
e
T
;
v
h
)
,
1
Δ
t
(
e
T
n
+
1
-
e
T
n
,
S
h
)
+
k
(
∇
e
T
n
+
1
/
2
,
∇
S
h
)
+
μ
h
(
∇
(
e
T
n
+
1
-
e
T
n
)
,
∇
S
h
)
=
E
N
n
(
u
,
T
;
S
h
)
,
with
(59)
M
O
n
(
u
,
p
,
T
,
e
T
;
v
h
)
=
(
u
(
t
n
+
1
)
-
u
(
t
n
)
Δ
t
-
u
t
(
t
n
+
1
/
2
)
,
v
h
)
+
Pr
(
∇
(
u
(
t
n
+
1
)
+
u
(
t
n
)
2
-
u
(
t
n
+
1
/
2
)
)
,
∇
v
h
)
+
μ
h
(
∇
(
u
(
t
n
+
1
)
-
u
(
t
n
)
)
,
∇
v
h
)
+
(
p
(
t
n
+
1
/
2
)
-
p
(
t
n
+
1
)
+
p
(
t
n
)
2
,
∇
·
v
h
)
+
Pr
Ra
(
ζ
e
T
n
+
1
/
2
,
v
h
)
-
Pr
Ra
(
ζ
T
(
t
n
+
1
)
+
T
(
t
n
)
2
,
v
h
)
+
Pr
Ra
(
ζ
T
(
t
n
+
1
/
2
)
,
v
h
)
-
c
(
u
(
t
n
+
1
/
2
)
,
u
(
t
n
+
1
/
2
)
,
v
h
)
+
c
(
E
[
u
h
n
,
u
h
n
-
1
]
,
u
h
n
+
1
+
u
h
n
2
,
v
h
)
,
E
N
n
(
u
,
T
;
S
h
)
=
(
T
(
t
n
+
1
)
-
T
(
t
n
)
Δ
t
-
T
t
(
t
n
+
1
/
2
)
,
S
h
)
+
k
(
∇
(
T
(
t
n
+
1
)
+
T
(
t
n
)
2
-
T
(
t
n
+
1
/
2
)
)
,
∇
S
h
)
+
μ
h
(
∇
(
T
(
t
n
+
1
)
-
T
(
t
n
)
)
,
∇
S
h
)
-
c
-
(
u
(
t
n
+
1
/
2
)
,
T
(
t
n
+
1
/
2
)
,
S
h
)
+
c
-
(
E
[
u
h
n
,
u
h
n
-
1
]
,
T
h
n
+
1
+
T
h
n
2
,
S
h
)
.
Using (31) and taking v
h
=
φ
u
n
+
1
/
2, S
h
=
φ
T
n
+
1
/
2 give
(60)
1
2
Δ
t
(
∥
φ
u
n
+
1
∥
2
-
∥
φ
u
n
∥
2
)
+
Pr
∥
∇
φ
u
n
+
1
/
2
∥
2
+
μ
h
2
(
∥
∇
φ
u
n
+
1
∥
2
-
∥
∇
φ
u
n
∥
2
)
=
(
η
u
n
+
1
-
η
u
n
Δ
t
,
φ
u
n
+
1
/
2
)
+
Pr
(
∇
η
u
n
+
1
/
2
,
∇
φ
u
n
+
1
/
2
)
-
(
e
p
n
+
1
/
2
,
∇
·
φ
u
n
+
1
/
2
)
+
μ
h
Δ
t
(
∇
(
η
u
n
+
1
-
η
u
n
Δ
t
)
,
∇
φ
u
n
+
1
/
2
)
-
M
O
n
(
u
,
p
,
T
,
e
T
;
φ
u
n
+
1
/
2
)
,
(61)
1
2
Δ
t
(
∥
φ
T
n
+
1
∥
2
-
∥
φ
T
n
∥
2
)
+
k
∥
∇
φ
T
n
+
1
/
2
∥
2
+
μ
h
2
(
∥
∇
φ
T
n
+
1
∥
2
-
∥
∇
φ
T
n
∥
2
)
=
(
η
T
n
+
1
-
η
T
n
Δ
t
,
φ
T
n
+
1
/
2
)
+
k
(
∇
η
T
n
+
1
/
2
,
∇
φ
T
n
+
1
/
2
)
+
μ
h
Δ
t
(
∇
(
η
T
n
+
1
-
η
T
n
Δ
t
)
,
∇
φ
T
n
+
1
/
2
)
-
E
N
n
(
u
,
T
;
φ
T
n
+
1
/
2
)
.
Step
2.2. Derive the error estimate of the momentum equation (60). Basing on Definition 2, we obtain
(62)
Pr
(
∇
η
u
n
+
1
/
2
,
∇
φ
u
n
+
1
/
2
)
-
(
e
p
n
+
1
/
2
,
∇
·
φ
u
n
+
1
/
2
)
=
0
,
k
(
∇
η
T
n
+
1
/
2
,
∇
φ
T
n
+
1
/
2
)
=
0
.
Similarly, for the linear terms on the RHS of (60), we have
(63)
|
(
η
u
n
+
1
-
η
u
n
Δ
t
,
φ
u
n
+
1
/
2
)
|
⩽
Pr
24
∥
∇
φ
u
n
+
1
/
2
∥
2
+
6
C
P
2
Pr
∥
η
u
n
+
1
-
η
u
n
Δ
t
∥
2
⩽
Pr
24
∥
∇
φ
u
n
+
1
/
2
∥
2
+
C
h
2
m
+
2
Pr
Δ
t
∥
u
t
∥
L
2
(
t
n
,
t
n
+
1
;
H
m
+
1
(
Ω
)
)
2
,
|
μ
h
Δ
t
(
∇
(
η
u
n
+
1
-
η
u
n
Δ
t
)
,
∇
φ
u
n
+
1
/
2
)
|
⩽
Pr
24
∥
∇
φ
u
n
+
1
/
2
∥
2
+
C
μ
2
Δ
t
h
2
m
+
2
Pr
∥
u
t
∥
L
2
(
t
n
,
t
n
+
1
;
H
m
+
1
(
Ω
)
)
2
.
For the linear terms of |
M
O
n
(
u
,
p
,
T
,
e
T
;
φ
u
n
+
1
/
2
)
| in (60), we obtain
(64)
|
(
u
(
t
n
+
1
)
-
u
(
t
n
)
Δ
t
-
u
t
(
t
n
+
1
/
2
)
,
φ
u
n
+
1
/
2
)
|
⩽
Pr
24
∥
∇
φ
u
n
+
1
/
2
∥
2
+
C
Δ
t
4
Pr
max
⁡
t
∈
[
t
n
,
t
n
+
1
]
∥
u
t
t
t
(
t
)
∥
2
,
Pr
|
(
∇
(
u
(
t
n
+
1
)
+
u
(
t
n
)
2
-
u
(
t
n
+
1
/
2
)
)
,
∇
φ
u
n
+
1
/
2
)
|
⩽
Pr
24
∥
∇
φ
u
n
+
1
/
2
∥
2
+
C
Pr
Δ
t
4
max
⁡
t
∈
[
t
n
,
t
n
+
1
]
∥
∇
u
t
t
(
t
)
∥
2
,
μ
h
Δ
t
|
(
∇
(
u
(
t
n
+
1
)
-
u
(
t
n
)
Δ
t
)
,
∇
φ
u
n
+
1
/
2
)
|
⩽
Pr
24
∥
∇
φ
u
n
+
1
/
2
∥
2
+
C
μ
2
h
2
Δ
t
2
Pr
max
⁡
t
∈
[
t
n
,
t
n
+
1
]
∥
∇
u
t
(
t
)
∥
2
,
|
(
p
(
t
n
+
1
)
+
p
(
t
n
)
2
-
p
(
t
n
+
1
/
2
)
,
∇
·
φ
u
n
+
1
/
2
)
|
⩽
Pr
24
∥
∇
φ
u
n
+
1
/
2
∥
2
+
C
Δ
t
4
Pr
max
⁡
t
∈
[
t
n
,
t
n
+
1
]
∥
p
t
t
(
t
)
∥
2
,
Pr
Ra
|
(
ζ
T
(
t
n
+
1
)
+
T
(
t
n
)
2
-
ζ
T
(
t
n
+
1
/
2
)
,
φ
u
n
+
1
/
2
)
|
⩽
Pr
24
∥
∇
φ
u
n
+
1
/
2
∥
2
+
C
Pr
R
a
2
Δ
t
4
max
⁡
t
∈
[
t
n
,
t
n
+
1
]
∥
T
t
t
(
t
)
∥
2
,
Pr
Ra
|
(
ζ
e
T
n
+
1
/
2
,
φ
u
n
+
1
/
2
)
|
=
Pr
Ra
|
(
ζ
η
T
n
+
1
/
2
-
ζ
φ
T
n
+
1
/
2
,
φ
u
n
+
1
/
2
)
|
⩽
Pr
Ra
(
C
P
∥
φ
T
n
+
1
/
2
∥
∥
∇
φ
u
n
+
1
/
2
∥
v
v
v
v
v
w
v
v
v
v
+
C
P
∥
η
T
n
+
1
/
2
∥
∥
∇
φ
u
n
+
1
/
2
∥
)
⩽
Pr
12
∥
∇
φ
u
n
+
1
/
2
∥
2
+
6
Pr
R
a
2
C
P
2
∥
φ
T
n
+
1
/
2
∥
2
+
6
Pr
R
a
2
C
P
2
∥
η
T
n
+
1
/
2
∥
2
.
For the nonlinear terms of |
M
O
n
(
u
,
p
,
T
,
e
T
;
φ
u
n
+
1
/
2
)
| in (60), it follows that
(65)
-
c
(
u
(
t
n
+
1
/
2
)
,
u
(
t
n
+
1
/
2
)
,
φ
u
n
+
1
/
2
)
+
c
(
E
[
u
h
n
,
u
h
n
-
1
]
,
u
h
n
+
1
+
u
h
n
2
,
φ
u
n
+
1
/
2
)
=
c
(
u
(
t
n
+
1
/
2
)
,
u
(
t
n
+
1
)
+
u
(
t
n
)
2
-
u
(
t
n
+
1
/
2
)
,
0000000
φ
u
n
+
1
/
2
u
(
t
n
+
1
)
+
u
(
t
n
)
2
)
+
c
(
u
(
t
n
+
1
)
+
u
(
t
n
)
2
E
[
u
(
t
n
)
,
u
(
t
n
-
1
)
]
00000000
-
u
(
t
n
+
1
/
2
)
,
u
(
t
n
+
1
)
+
u
(
t
n
)
2
,
φ
u
n
+
1
/
2
)
-
c
(
E
[
η
u
n
,
η
u
n
-
1
]
,
u
(
t
n
+
1
)
+
u
(
t
n
)
2
,
φ
u
n
+
1
/
2
)
+
c
(
E
[
φ
u
n
,
φ
u
n
-
1
]
,
u
(
t
n
+
1
)
+
u
(
t
n
)
2
,
φ
u
n
+
1
/
2
)
-
c
(
E
[
u
h
n
,
u
h
n
-
1
]
,
η
u
n
+
1
/
2
,
φ
u
n
+
1
/
2
)
=
I
1
+
I
2
+
⋯
+
I
5
,
where we have used Lemma 1 and rearranged some terms. Using definition (15) of the operator E
[
·
,
·
] and the regularity assumptions (32) yields
(66)
∥
∇
E
[
u
(
t
n
)
,
u
(
t
n
-
1
)
]
∥
⩽
C
,
∥
∇
E
[
η
u
n
,
η
u
n
-
1
]
∥
⩽
3
2
∥
∇
η
u
n
∥
+
1
2
∥
∇
η
u
n
-
1
∥
.
We now estimate each nonlinear terms of (65). For the terms I
1
,
…
,
I
4, by making use of Lemma 1, Taylor expansion on t and Young’s inequality together with the regularity assumptions (32), we infer that
(67)
|
I
1
|
+
|
I
2
|
⩽
C
∥
∇
(
u
(
t
n
+
1
)
+
u
(
t
n
)
2
-
u
(
t
n
+
1
/
2
)
)
∥
×
∥
∇
u
(
t
n
+
1
/
2
)
∥
∥
∇
φ
u
n
+
1
/
2
∥
+
C
∥
∇
(
3
2
u
(
t
n
)
-
1
2
u
(
t
n
)
-
u
(
t
n
+
1
/
2
)
)
∥
×
∥
∇
(
u
(
t
n
+
1
)
+
u
(
t
n
)
2
)
∥
∥
∇
φ
u
n
+
1
/
2
∥
⩽
Pr
24
∥
∇
φ
u
n
+
1
/
2
∥
2
+
C
Δ
t
4
Pr
max
⁡
t
∈
[
t
n
,
t
n
+
1
]
∥
∇
u
t
t
(
t
)
∥
2
,
|
I
3
|
⩽
C
∥
∇
E
[
η
u
n
,
η
u
n
-
1
]
∥
∥
∇
φ
u
n
+
1
/
2
∥
⩽
Pr
24
∥
∇
φ
u
n
+
1
/
2
∥
2
+
C
P
r
-
1
(
∥
∇
η
u
n
∥
2
+
∥
∇
η
u
n
-
1
∥
2
)
,
|
I
4
|
⩽
C
∥
E
[
φ
u
n
,
φ
u
n
-
1
]
∥
∥
∇
φ
u
n
+
1
/
2
∥
⩽
Pr
48
∥
∇
φ
u
n
+
1
/
2
∥
2
+
C
P
r
-
1
(
∥
φ
u
n
∥
2
+
∥
φ
u
n
-
1
∥
2
)
.
Adding and subtracting c
(
E
[
u
(
t
n
)
,
u
(
t
n
-
1
)
]
,
η
u
n
+
1
/
2
,
φ
u
n
+
1
/
2
) to I
5, using Lemma 1, the inverse, and Young’s inequalities, we have (for detail, see [25])
(68)
|
I
5
|
⩽
Pr
48
∥
∇
φ
u
n
+
1
/
2
∥
2
+
C
P
r
-
1
(
1
+
∥
∇
η
u
n
∥
2
+
∥
∇
η
u
n
-
1
∥
2
)
∥
∇
η
u
n
+
1
/
2
∥
2
+
C
h
-
d
/
2
∥
∇
η
u
n
+
1
/
2
∥
(
∥
φ
u
n
-
1
∥
2
+
∥
φ
u
n
∥
2
+
∥
φ
u
n
+
1
∥
2
)
.
Combining these inequalities with (60), we arrive at
(69)
1
2
Δ
t
(
∥
φ
u
n
+
1
∥
2
-
∥
φ
u
n
∥
2
)
+
Pr
2
∥
∇
φ
u
n
+
1
/
2
∥
2
+
μ
h
2
(
∥
∇
φ
u
n
+
1
∥
2
-
∥
∇
φ
u
n
∥
2
)
⩽
C
h
2
m
+
2
Pr
(
1
Δ
t
+
μ
2
Δ
t
)
∥
u
t
∥
L
2
(
t
n
,
t
n
+
1
;
H
m
+
1
(
Ω
)
)
2
+
C
Δ
t
4
Pr
[
max
⁡
t
∈
[
t
n
,
t
n
+
1
]
∥
u
t
t
t
(
t
)
∥
2
000000000000
+
(
1
+
P
r
2
)
max
⁡
t
∈
[
t
n
,
t
n
+
1
]
∥
∇
u
t
t
(
t
)
∥
2
000000000000
+
P
r
2
R
a
2
max
⁡
t
∈
[
t
n
,
t
n
+
1
]
∥
T
t
t
(
t
)
∥
2
000000000000
+
max
⁡
t
∈
[
t
n
,
t
n
+
1
]
∥
p
t
t
(
t
)
∥
2
]
+
C
μ
2
h
2
Δ
t
2
Pr
max
⁡
t
∈
[
t
n
,
t
n
+
1
]
∥
∇
u
t
(
t
)
∥
2
+
C
P
r
-
1
(
∥
∇
η
u
n
∥
2
+
∥
∇
η
u
n
-
1
∥
2
)
+
C
P
r
-
1
(
∥
∇
φ
u
n
∥
2
+
∥
∇
φ
u
n
-
1
∥
2
)
+
6
Pr
R
a
2
C
P
2
(
∥
η
T
n
+
1
/
2
∥
2
+
∥
φ
T
n
+
1
/
2
∥
2
)
+
C
P
r
-
1
(
1
+
∥
∇
η
u
n
∥
2
+
∥
∇
η
u
n
-
1
∥
2
)
∥
∇
η
u
n
+
1
/
2
∥
2
+
C
h
-
d
/
2
∥
∇
η
u
n
+
1
/
2
∥
(
∥
φ
u
n
-
1
∥
2
+
∥
φ
u
n
∥
2
+
∥
φ
u
n
+
1
∥
2
)
.
Step
2.3. Establish the error estimate of the energy error equation (61). For the linear terms on the RHS of (61), we yield that
(70)
|
(
η
T
n
+
1
-
η
T
n
Δ
t
,
φ
T
n
+
1
/
2
)
|
⩽
k
min
⁡
18
∥
∇
φ
T
n
+
1
/
2
∥
2
+
9
C
P
2
k
min
⁡
∥
η
T
n
+
1
-
η
T
n
Δ
t
∥
2
⩽
k
min
⁡
18
∥
∇
φ
T
n
+
1
/
2
∥
2
+
C
h
2
m
+
2
k
min
⁡
Δ
t
∥
T
t
∥
L
2
(
t
n
,
t
n
+
1
;
H
m
+
1
(
Ω
)
)
2
,
|
μ
h
Δ
t
(
∇
(
η
T
n
+
1
-
η
T
n
Δ
t
)
,
∇
φ
T
n
+
1
/
2
)
|
⩽
k
min
⁡
18
∥
∇
φ
T
n
+
1
/
2
∥
2
+
C
μ
2
Δ
t
h
2
m
+
2
k
min
⁡
∥
T
t
∥
L
2
(
t
n
,
t
n
+
1
;
H
m
+
1
(
Ω
)
)
2
.
For the linear terms of|
E
N
n
(
u
,
T
;
φ
T
n
+
1
/
2
)
|, we obtain
(71)
|
(
T
(
t
n
+
1
)
-
T
(
t
n
)
Δ
t
-
T
t
(
t
n
+
1
/
2
)
,
φ
T
n
+
1
/
2
)
|
⩽
k
min
⁡
18
∥
∇
φ
T
n
+
1
/
2
∥
2
+
C
Δ
t
4
k
min
⁡
max
⁡
t
∈
[
t
n
,
t
n
+
1
]
∥
T
t
t
t
(
t
)
∥
2
,
k
|
(
∇
(
T
(
t
n
+
1
)
+
T
(
t
n
)
2
-
T
(
t
n
+
1
/
2
)
)
,
∇
φ
T
n
+
1
/
2
)
|
⩽
k
min
⁡
18
∥
∇
φ
T
n
+
1
/
2
∥
2
+
C
Δ
t
4
k
max
⁡
2
k
min
⁡
max
⁡
t
∈
[
t
n
,
t
n
+
1
]
∥
∇
T
t
t
(
t
)
∥
2
,
μ
h
Δ
t
|
(
∇
(
T
(
t
n
+
1
)
-
T
(
t
n
)
Δ
t
)
,
∇
φ
T
n
+
1
/
2
)
|
⩽
k
min
⁡
18
∥
∇
φ
T
n
+
1
/
2
∥
2
+
C
μ
2
h
2
Δ
t
2
k
min
⁡
max
⁡
t
∈
[
t
n
,
t
n
+
1
]
∥
∇
T
t
(
t
)
∥
2
.
By adding and subtracting c
(
E
[
u
h
n
,
u
h
n
-
1
]
-
E
[
u
(
t
n
)
,
u
(
t
n
-
1
)
]
,
(
T
(
t
n
+
1
)
+
T
(
t
n
)
)
/
2
,
φ
T
n
+
1
/
2
) to the nonlinear terms of |
E
N
n
(
u
,
T
,
φ
T
n
+
1
/
2
)
| in (61), we can obtain(72)
-
c
-
(
u
(
t
n
+
1
/
2
)
,
T
(
t
n
+
1
/
2
)
,
φ
T
n
+
1
/
2
)
+
c
-
(
E
[
u
h
n
,
u
h
n
-
1
]
,
T
h
n
+
1
+
T
h
n
2
,
φ
T
n
+
1
/
2
)
=
c
-
(
u
(
t
n
+
1
/
2
)
,
T
(
t
n
+
1
)
+
T
(
t
n
)
2
0000000
-
T
(
t
n
+
1
/
2
)
,
φ
T
n
+
1
/
2
u
(
t
n
+
1
)
+
u
(
t
n
)
2
)
+
c
-
(
u
(
t
n
+
1
)
+
u
(
t
n
)
2
E
[
u
(
t
n
)
,
u
(
t
n
-
1
)
]
0000000
-
u
(
t
n
+
1
/
2
)
,
T
(
t
n
+
1
)
+
T
(
t
n
)
2
,
φ
T
n
+
1
/
2
)
-
c
-
(
E
[
η
u
n
,
η
u
n
-
1
]
,
T
(
t
n
+
1
)
+
T
(
t
n
)
2
,
φ
T
n
+
1
/
2
)
+
c
-
(
E
[
φ
u
n
,
φ
u
n
-
1
]
,
T
(
t
n
+
1
)
+
T
(
t
n
)
2
,
φ
T
n
+
1
/
2
)
-
c
-
(
E
[
u
h
n
,
u
h
n
-
1
]
,
η
T
n
+
1
/
2
,
φ
T
n
+
1
/
2
)
=
I
I
1
+
I
I
2
+
⋯
+
I
I
5
.
Same as Step 2.2, basing on Lemma 1, Taylor expansion on t, Young’s inequality, the inverse inequality, the regularity assumptions (32), and the estimate (66) of the operator E
[
·
,
·
], we can derive that
(73)
|
I
I
1
|
+
|
I
I
2
|
⩽
C
∥
∇
(
T
(
t
n
+
1
)
+
T
(
t
n
)
2
-
T
(
t
n
+
1
/
2
)
)
∥
×
∥
∇
u
(
t
n
+
1
/
2
)
∥
∥
∇
φ
T
n
+
1
/
2
∥
+
C
∥
∇
(
3
2
u
(
t
n
)
-
1
2
u
(
t
n
)
-
u
(
t
n
+
1
/
2
)
)
∥
×
∥
∇
(
T
(
t
n
+
1
)
+
T
(
t
n
)
2
)
∥
∥
∇
φ
T
n
+
1
/
2
∥
⩽
k
min
⁡
18
∥
∇
φ
T
n
+
1
/
2
∥
2
+
C
Δ
t
4
k
min
⁡
(
max
⁡
t
∈
[
t
n
,
t
n
+
1
]
∥
∇
u
t
t
(
t
)
∥
2
+
max
⁡
t
∈
[
t
n
,
t
n
+
1
]
∥
∇
T
t
t
(
t
)
∥
2
)
,
|
I
I
3
|
⩽
C
∥
∇
E
[
η
u
n
,
η
u
n
-
1
]
∥
∥
∇
φ
T
n
+
1
/
2
∥
⩽
k
min
⁡
18
∥
∇
φ
T
n
+
1
/
2
∥
2
+
C
k
min
⁡
-
1
(
∥
∇
η
u
n
∥
2
+
∥
∇
η
u
n
-
1
∥
2
)
,
|
I
I
4
|
⩽
C
∥
E
[
φ
u
n
,
φ
u
n
-
1
]
∥
∥
∇
φ
T
n
+
1
/
2
∥
⩽
k
min
⁡
18
∥
∇
φ
T
n
+
1
/
2
∥
2
+
C
k
min
⁡
-
1
(
∥
φ
u
n
∥
2
+
∥
φ
u
n
-
1
∥
2
)
,
|
I
I
5
|
⩽
k
min
⁡
18
∥
∇
φ
T
n
+
1
/
2
∥
2
+
C
k
min
⁡
-
1
∥
∇
η
T
n
+
1
/
2
∥
2
(
1
+
∥
∇
η
u
n
∥
2
+
∥
∇
η
u
n
-
1
∥
2
)
+
C
h
-
d
/
2
∥
∇
η
T
n
+
1
/
2
∥
×
(
∥
φ
u
n
-
1
∥
2
+
∥
φ
u
n
∥
2
+
∥
φ
T
n
∥
2
+
∥
φ
T
n
+
1
∥
2
)
.
Substituting the above inequalities of Step 2.3 into (61), we obtain
(74)
1
2
Δ
t
(
∥
φ
T
n
+
1
∥
2
-
∥
φ
T
n
∥
2
)
+
k
min
⁡
2
∥
∇
φ
T
n
+
1
/
2
∥
2
vvvv
+
μ
h
2
(
∥
∇
φ
T
n
+
1
∥
2
-
∥
∇
φ
T
n
∥
2
)
⩽
C
h
2
m
+
2
k
min
⁡
(
1
Δ
t
+
μ
2
Δ
t
)
∥
T
t
∥
L
2
(
t
n
,
t
n
+
1
;
H
m
+
1
(
Ω
)
)
2
vvvv
+
C
Δ
t
4
k
min
⁡
[
max
⁡
t
∈
[
t
n
,
t
n
+
1
]
∥
T
t
t
t
(
t
)
∥
2
000000000000
+
(
1
+
k
max
⁡
2
)
max
⁡
t
∈
[
t
n
,
t
n
+
1
]
∥
∇
T
t
t
(
t
)
∥
2
000000000000
+
max
⁡
t
∈
[
t
n
,
t
n
+
1
]
∥
∇
u
t
t
(
t
)
∥
2
]
vvvv
+
C
μ
2
h
2
Δ
t
2
k
min
⁡
max
⁡
t
∈
[
t
n
,
t
n
+
1
]
∥
∇
T
t
(
t
)
∥
2
vvvv
+
C
k
min
⁡
-
1
(
∥
∇
η
u
n
∥
2
+
∥
∇
η
u
n
-
1
∥
2
)
vvvv
+
C
k
min
⁡
-
1
(
∥
φ
u
n
∥
2
+
∥
φ
u
n
-
1
∥
2
)
vvvv
+
C
k
min
⁡
-
1
(
1
+
∥
∇
η
u
n
∥
2
+
∥
∇
η
u
n
-
1
∥
2
)
∥
∇
η
T
n
+
(
1
/
2
)
∥
2
vvvv
+
C
h
-
d
/
2
∥
∇
η
T
n
+
1
/
2
∥
(
∥
φ
u
n
-
1
∥
2
+
∥
φ
u
n
∥
2
000000000000000000000
+
∥
φ
T
n
∥
2
+
∥
φ
T
n
+
1
∥
2
)
.
Step
2.4. Combining (69) with (74) and rearranging some terms yield
(75)
1
2
Δ
t
(
∥
φ
u
n
+
1
∥
2
-
∥
φ
u
n
∥
2
)
+
Pr
2
∥
∇
φ
u
n
+
1
/
2
∥
2
vvvv
+
μ
h
2
(
∥
∇
φ
u
n
+
1
∥
2
-
∥
∇
φ
u
n
∥
2
)
+
1
2
Δ
t
(
∥
φ
T
n
+
1
∥
2
-
∥
φ
T
n
∥
2
)
vvvv
+
k
min
⁡
2
∥
∇
φ
T
n
+
1
/
2
∥
2
+
μ
h
2
(
∥
∇
φ
T
n
+
1
∥
2
-
∥
∇
φ
T
n
∥
2
)
⩽
C
h
2
m
+
2
(
1
Δ
t
+
μ
2
Δ
t
)
vvvv
×
[
∥
u
t
∥
L
2
(
t
n
,
t
n
+
1
;
H
m
+
1
(
Ω
)
)
2
Pr
+
∥
T
t
∥
L
2
(
t
n
,
t
n
+
1
;
H
m
+
1
(
Ω
)
)
2
k
min
⁡
]
vvvv
+
C
Δ
t
4
[
P
r
-
1
max
⁡
t
∈
[
t
n
,
t
n
+
1
]
∥
u
t
t
t
(
t
)
∥
2
000000000000
+
(
Pr
+
P
r
-
1
+
k
min
⁡
-
1
)
max
⁡
t
∈
[
t
n
,
t
n
+
1
]
∥
∇
u
t
t
(
t
)
∥
2
000000000000
+
P
r
-
1
max
⁡
t
∈
[
t
n
,
t
n
+
1
]
∥
p
t
t
(
t
)
∥
2
000000000000
+
P
r
2
R
a
2
max
⁡
t
∈
[
t
n
,
t
n
+
1
]
∥
T
t
t
(
t
)
∥
2
000000000000
+
k
min
⁡
-
1
max
⁡
t
∈
[
t
n
,
t
n
+
1
]
∥
T
t
t
t
(
t
)
∥
2
000000000000
+
k
min
⁡
-
1
(
1
+
k
max
⁡
2
)
max
⁡
t
∈
[
t
n
,
t
n
+
1
]
∥
∇
T
t
t
(
t
)
∥
2
]
vvvv
+
C
μ
2
h
2
Δ
t
2
[
P
r
-
1
max
⁡
t
∈
[
t
n
,
t
n
+
1
]
∥
∇
u
t
(
t
)
∥
2
000000000000000
+
k
min
⁡
-
1
max
⁡
t
∈
[
t
n
,
t
n
+
1
]
∥
∇
T
t
(
t
)
∥
2
]
vvvv
+
C
(
P
r
-
1
+
k
min
⁡
-
1
)
vvvv
×
(
∥
∇
η
u
n
∥
2
+
∥
∇
η
u
n
-
1
∥
2
+
∥
φ
u
n
∥
2
+
∥
φ
u
n
-
1
∥
2
)
vvvv
+
C
(
1
+
∥
∇
η
u
n
∥
2
+
∥
∇
η
u
n
-
1
∥
2
)
vvvv
×
(
P
r
-
1
∥
∇
η
u
n
+
1
/
2
∥
2
+
k
min
⁡
-
1
∥
∇
η
T
n
+
1
/
2
∥
2
)
vvvv
+
C
h
-
d
/
2
∥
∇
η
u
n
+
1
/
2
∥
vvvv
×
(
∥
φ
u
n
-
1
∥
2
+
∥
φ
u
n
∥
2
+
∥
φ
u
n
+
1
∥
2
)
vvvv
+
C
h
-
d
/
2
∥
∇
η
T
n
+
1
/
2
∥
vvvv
×
(
∥
φ
u
n
-
1
∥
2
+
∥
φ
u
n
∥
2
+
∥
φ
T
n
∥
2
+
∥
φ
T
n
+
1
∥
2
)
vvvv
+
6
Pr
R
a
2
C
P
2
(
∥
η
T
n
+
1
/
2
∥
2
+
∥
φ
T
n
+
1
/
2
∥
2
)
.
Multiplying (75) by 2
Δ
t and summing it over n from 1 to l, using the approximation property (9), and substituting the error estimate (55) of the first time level into it, we can obtain
(76)
∥
φ
u
l
+
1
∥
2
+
Δ
t
∑
n
=
0
l
Pr
∥
∇
φ
u
l
+
1
/
2
∥
2
+
μ
h
Δ
t
∥
∇
φ
u
l
+
1
∥
2
+
∥
φ
T
l
+
1
∥
2
+
Δ
t
∑
n
=
0
l
k
min
⁡
∥
∇
φ
T
l
+
1
/
2
∥
2
+
μ
h
Δ
t
∥
∇
φ
T
l
+
1
∥
2
⩽
C
h
2
m
+
2
(
1
+
μ
2
Δ
t
2
)
×
[
∥
u
t
∥
L
2
(
0
,
t
*
;
H
m
+
1
(
Ω
)
)
2
Pr
+
∥
T
t
∥
L
2
(
0
,
t
*
;
H
m
+
1
(
Ω
)
)
2
k
min
⁡
]
+
C
Δ
t
5
[
P
r
-
1
∥
u
t
t
t
(
t
)
∥
L
∞
(
0
,
t
*
;
L
2
(
Ω
)
)
2
00000000000
+
P
r
-
1
∥
p
t
t
(
t
)
∥
L
∞
(
0
,
t
*
;
L
2
(
Ω
)
)
2
00000000000
+
(
Pr
+
P
r
-
1
+
k
min
⁡
-
1
)
∥
∇
u
t
t
(
t
)
∥
L
∞
(
0
,
t
*
;
L
2
(
Ω
)
)
2
00000000000
+
P
r
2
R
a
2
∥
T
t
t
(
t
)
∥
L
∞
(
0
,
t
*
;
L
2
(
Ω
)
)
2
00000000000
+
k
min
⁡
-
1
∥
T
t
t
t
(
t
)
∥
L
∞
(
0
,
t
*
;
L
2
(
Ω
)
)
2
00000000000
+
k
min
⁡
-
1
(
1
+
k
max
⁡
2
)
∥
∇
T
t
t
(
t
)
∥
L
∞
(
0
,
t
*
;
L
2
(
Ω
)
)
2
]
+
C
μ
2
h
2
Δ
t
3
[
P
r
-
1
∥
∇
u
t
(
t
)
∥
L
∞
(
0
,
t
*
;
L
2
(
Ω
)
)
2
000000000000000
+
k
min
⁡
-
1
∥
∇
T
t
(
t
)
∥
L
∞
(
0
,
t
*
;
L
2
(
Ω
)
)
2
]
+
C
(
P
r
-
1
+
k
min
⁡
-
1
)
Δ
t
h
2
m
∥
u
∥
L
2
(
0
,
t
*
;
H
m
+
1
(
Ω
)
)
2
+
C
k
min
⁡
-
1
Δ
t
h
2
m
∥
T
∥
L
2
(
0
,
t
*
;
H
m
+
1
(
Ω
)
)
2
+
C
(
P
r
-
1
+
k
min
⁡
-
1
)
Δ
t
∑
n
=
0
l
∥
φ
u
n
∥
2
+
C
h
m
-
(
d
/
2
)
∥
u
∥
L
∞
(
0
,
t
*
;
H
m
+
1
(
Ω
)
)
Δ
t
∥
φ
u
l
+
1
∥
2
+
C
h
m
-
(
d
/
2
)
(
∥
u
∥
L
∞
(
0
,
t
*
;
H
m
+
1
(
Ω
)
)
00000000000000
+
∥
T
∥
L
∞
(
0
,
t
*
;
H
m
+
1
(
Ω
)
)
)
Δ
t
∑
n
=
0
l
∥
φ
u
n
∥
2
+
(
C
h
m
-
(
d
/
2
)
∥
T
∥
L
∞
(
0
,
t
*
;
H
m
+
1
(
Ω
)
)
000000
+
6
Pr
R
a
2
C
P
2
)
Δ
t
∑
n
=
0
l
∥
φ
T
l
∥
2
+
(
C
h
m
-
(
d
/
2
)
∥
T
∥
L
∞
(
0
,
t
*
;
H
m
+
1
(
Ω
)
)
Δ
t
000000
+
6
Pr
R
a
2
C
P
2
Δ
t
)
∥
φ
T
l
+
1
∥
2
+
C
Pr
R
a
2
C
P
2
h
2
m
+
2
Δ
t
∥
T
∥
L
2
(
0
,
t
*
;
H
m
+
1
(
Ω
)
)
2
.
Applying condition (33) and the regularity assumptions (32) to (76) yields
(77)
1
2
∥
φ
u
l
+
1
∥
2
+
Δ
t
∑
n
=
0
l
Pr
∥
∇
φ
u
l
+
1
/
2
∥
2
+
μ
h
Δ
t
∥
∇
φ
u
l
+
1
∥
2
+
1
2
∥
φ
T
l
+
1
∥
2
+
Δ
t
∑
n
=
0
l
k
min
⁡
∥
∇
φ
T
l
+
(
1
/
2
)
∥
2
+
μ
h
Δ
t
∥
∇
φ
T
l
+
1
∥
2
⩽
C
(
Pr
+
P
r
-
1
+
k
min
⁡
-
1
+
k
min
⁡
-
1
k
max
⁡
2
+
Pr
R
a
2
C
P
2
)
×
(
h
2
m
+
μ
2
h
2
Δ
t
2
+
Δ
t
4
)
+
C
(
h
m
-
(
d
/
2
)
∥
T
∥
L
∞
(
0
,
t
*
;
H
m
+
1
(
Ω
)
)
00000000
+
6
Pr
R
a
2
C
P
2
h
m
-
(
d
/
2
)
∥
T
∥
L
∞
(
0
,
t
*
;
H
m
+
1
(
Ω
)
)
)
Δ
t
∑
n
=
0
l
∥
φ
T
n
∥
2
+
C
(
P
r
-
1
+
k
min
⁡
-
1
000000000
+
h
m
-
(
d
/
2
)
(
∥
u
∥
L
∞
(
0
,
t
*
;
H
m
+
1
(
Ω
)
)
000000000000000000
+
∥
T
∥
L
∞
(
0
,
t
*
;
H
m
+
1
(
Ω
)
)
P
r
-
1
)
)
Δ
t
∑
n
=
0
l
∥
φ
u
n
∥
2
.
Making use of the Gronwall lemma [22], for any 0
⩽
l
⩽
N
-
1, there exists a constant C
^
=
C
^
(
P
r
,
k
,
R
a
,
Ω
,
m
,
u
,
p
,
T
) such that
(78)
1
2
∥
φ
u
l
+
1
∥
2
+
Δ
t
∑
n
=
0
l
P
r
∥
∇
φ
u
l
+
1
/
2
∥
2
+
μ
h
Δ
t
∥
∇
φ
u
l
+
1
∥
2
+
1
2
∥
φ
T
l
+
1
∥
2
+
Δ
t
∑
n
=
0
l
k
min
⁡
∥
∇
φ
T
l
+
(
1
/
2
)
∥
2
+
μ
h
Δ
t
∥
∇
φ
T
l
+
1
∥
2
⩽
C
^
(
h
2
m
+
μ
2
h
2
Δ
t
2
+
Δ
t
4
)
.
We complete the proof of Theorem 6 by using the triangle inequality.Corollary 8.
Under the conditions of Theorem6 and letting (
X
h
,
M
h
,
W
h
) be continuous (
P
2
,
P
1
,
P
2
) finite element spaces, then there exists a constant C
^
=
C
^
(
P
r
,
k
,
R
a
,
Ω
,
u
,
p
,
T
) such that, for any 0
⩽
l
⩽
N
-
1,
(79)
∥
u
(
t
l
+
1
)
-
u
h
l
+
1
∥
+
(
Δ
t
∑
n
=
0
l
P
r
∥
∇
(
(
(
u
(
t
n
+
1
)
-
u
h
n
+
1
)
+
(
u
(
t
n
)
-
u
h
n
)
)
000000000000000000
×
(
2
)
-
1
(
(
u
(
t
n
+
1
)
-
u
h
n
+
1
)
+
(
u
(
t
n
)
-
u
h
n
)
)
)
∥
2
∑
n
=
0
l
)
1
/
2
+
∥
T
(
t
l
+
1
)
-
T
h
l
+
1
∥
+
(
Δ
t
∑
n
=
0
l
k
min
⁡
∥
∇
(
(
(
T
(
t
n
+
1
)
-
T
h
n
+
1
)
+
(
T
(
t
n
)
-
T
h
n
)
)
000000000000000000
×
(
2
)
-
1
(
(
T
(
t
n
+
1
)
-
T
h
n
+
1
)
+
(
T
(
t
n
)
-
T
h
n
)
)
)
∥
2
∑
n
=
0
l
)
1
/
2
+
(
μ
h
Δ
t
)
1
/
2
∥
∇
(
u
(
t
l
+
1
)
-
u
h
l
+
1
)
∥
+
(
μ
h
Δ
t
)
1
/
2
∥
∇
(
T
(
t
l
+
1
)
-
T
h
l
+
1
)
∥
⩽
C
^
(
h
2
+
μ
h
Δ
t
+
Δ
t
2
)
.Remark 9.
Corollary8 shows that the optimal convergence rate of Algorithm 4 in space by using (
P
2
,
P
1
,
P
2
) elements is O
(
h
2
) for temperature and velocity in H
1 seminorm, with the choosing μ
=
O
(
1
) and Δ
t
=
O
(
h
). However, the error estimates in space are suboptimal for temperature and velocity in L
2 norm. We can find these facts from Tables 1–4 in the next section.Table 1
Errors and CPU costs of velocityu by using CNSLE and Δ
t
=
(
1
/
10
)
h.
1
/
h
∥
u
-
u
h
∥
L
∞
,
0
CR
∥
u
-
u
h
∥
L
0,1
CR
CPU
4
0.00166806
—
0.0149995
—
3.104
8
0.000194755
3.09844
0.00402588
1.89754
10.576
16
2.3595
e
-
005
3.04511
0.00103018
1.96641
52.494
32
2.93692
e
-
006
3.00611
0.000259323
1.99007
405.008
64
3.69237
e
-
007
2.99168
6.49529
e
-
005
1.99728
3274.88Table 2
Errors and CPU costs of temperatureT by using CNSLE and Δ
t
=
(
1
/
10
)
h.
1
/
h
∥
T
-
T
h
∥
L
∞
,
0
CR
∥
T
-
T
h
∥
L
0,1
CR
CPU
4
0.000922881
—
0.00869632
—
3.104
8
0.00011998
2.94335
0.00240038
1.85714
10.576
16
1.49476
e
-
005
3.00481
0.000615762
1.96282
52.494
32
1.86956
e
-
006
2.99914
0.000154998
1.99012
405.008
64
2.37001
e
-
007
2.97973
3.88187
e
-
005
1.99743
3274.88Table 3
Errors and CPU costs of velocityu by using CNS and Δ
t
=
(
1
/
10
)
h.
1
/
h
∥
u
-
u
h
∥
L
∞
,
0
CR
∥
u
-
u
h
∥
L
0,1
CR
CPU
4
0.00166784
—
0.0149995
—
8.174
8
0.000194742
3.09834
0.00402588
1.89754
19.609
16
2.35941
e
-
005
3.04507
0.00103018
1.96641
105.721
32
2.93687
e
-
006
3.00607
0.000259323
1.99007
752.208
64
3.69237
e
-
007
2.99166
6.49529
e
-
005
1.99728
5581.34Table 4
Errors and CPU costs of temperatureT by using CNS and Δ
t
=
(
1
/
10
)
h.
1
/
h
∥
T
-
T
h
∥
L
∞
,
0
CR
∥
T
-
T
h
∥
L
0,1
CR
CPU
4
0.000922881
—
0.00869632
—
8.174
8
0.00011998
2.94335
0.00240038
1.85714
19.609
16
1.49476
e
-
005
3.00481
0.000615762
1.96282
105.721
32
1.86956
e
-
006
2.99914
0.000154998
1.99012
752.208
64
2.37001
e
-
007
2.97973
3.88187
e
-
005
1.99743
5581.34


## 5. Numerical Experiments

We first verify the convergence rates and the effectiveness of our methods by an analytical solution. Then we test the stability of the method in caseRa is large in a squared cavity with the left wall heating. The code was implemented using the software package FreeFEM++ [37]. The pairs of continuous (
P
2
,
P
1
,
P
2
) elements are chosen for the FE spaces (
X
h
,
M
h
,
W
h
). All computations are carried out in the domain Ω
=
[
0,1
]
2. The uniform mesh is obtained by dividing Ω into squares and then drawing a diagonal in each square in the same direction. We set the parameter μ
=
O
(
1
).Example 1 (analytical solution).
As in [38–41], to obtain an analytical solution for the considered problem, a right-hand side function is added to the momentum equation of (1). The analytical solution is
(80)
u
1
(
x
,
y
)
=
10
x
2
(
x
-
1
)
2
y
(
y
-
1
)
(
2
y
-
1
)
cos
⁡
(
t
)
,
u
2
(
x
,
y
)
=
-
10
x
(
x
-
1
)
(
2
x
-
1
)
y
2
(
y
-
1
)
2
cos
⁡
(
t
)
,
p
(
x
,
y
)
=
10
(
2
x
-
1
)
(
2
y
-
1
)
cos
⁡
(
t
)
,
T
(
x
,
y
)
=
u
1
(
x
,
y
)
+
u
2
(
x
,
y
)
.
The initial velocity field and temperature field are equal to the analytical solution at time t
=
0. We choose the parameters that are k
=
1.0, Pr
=
1.0, and Ra
=
100, respectively. The mesh and time step sizes scalings are set that for a refinement, each of h and Δ
t get cut in half, where the final time is chosen as t
*
=
0.1 and Δ
t
=
(
1
/
10
)
h. We list the errors and the convergence rates (CR) for velocity u and temperature T in L
∞
(
0
,
t
*
;
L
2
(
Ω
)
) (denoting as L
∞
,
0) and L
2
(
0
,
t
*
;
H
1
(
Ω
)
) (denoting as L
0,1) norm in Tables 1–4, respectively. Tables 1 and 2 are the results of Crank-Nicolson linearized extrapolation stabilized (CNSLE) FE method. Tables 3 and 4 are the results of Crank-Nicolson stabilized (CNS) FE method, which the nonlinear system is solved by Newton iterative at each time step and the iterative tolerance is 1
0
-
8. We find that the errors listed in Tables 1 and 3 and in Tables 2 and 4 are similar, respectively, which means that the CNSLE method is comparable to the CNS method. From Tables 1–4, we find that a quadratic convergence rate for the computed velocity and temperature in H
1 seminorm is obtained, which test and verify the theoretical results. However in Tables 1–4, it is easy to see that a cubic convergence rate both for the computed velocity and temperature in L
2 norm is obtained, which indicate that our error estimate in L
2 norm is suboptimal. As expected, since CNSLE is the linearized version of the CNS method, which does not include any iterations in computing, CPU cost by CNSLE method is relatively less than by CNS method. The numerical results agree well with the theoretical predictions.Example 2 (natural convection in a squared cavity with the left wall heating).
Now we present natural convection in a squared cavity with the left wall heating. The boundary conditions are given in Figure1. We choose the initial conditions u
0
=
1, T
0
=
1, the parameters k
=
1.0, Pr
=
0.71, 1
0
4
⩽
Ra
⩽
1
0
7, the right hand functions γ
=
0, the time step Δ
t
=
0.1, and the uniform mesh size h
=
1
/
64. We performed the following study: calculating the problem from the rest until the approximation solution reaches a steady state. The criterion to stop this process is
(81)
max
⁡
{
∥
u
h
n
+
1
-
u
h
n
∥
L
2
(
Ω
)
∥
u
h
n
+
1
∥
L
2
(
Ω
)
,
∥
T
h
n
+
1
-
T
h
n
∥
L
2
(
Ω
)
∥
T
h
n
+
1
∥
L
2
(
Ω
)
}
⩽
1
0
-
5
,
where n
+
1, n denote t
n
+
1, t
n, respectively. We describe the final results of the problem at its steady state in Tables 5–7 and Figures 2–4, which are according to the results of [16, 17, 34–36]. The numbers in parenthesis of Tables 5–7 correspond to the used mesh.
Tables5 and 6 present the maximum vertical velocity at midheight and horizontal velocity at midwidth and compare the quantitative results with the available benchmark solutions [16, 34–36]. Figure 2 presents the vertical velocity distribution at the midheight and the horizontal velocity distribution at the midwidth for Ra
=
1
0
4, 1
0
5, 1
0
6, 1
0
7, respectively. We find that the differences in the profiles are getting larger along with the increase of Rayleigh numbers. It is seen that the agreement is excellent even at higher Rayleigh numbers.
The heat transfer coefficient in terms of the local Nusselt number is defined by(82)
Nu
local
(
x
,
y
)
=
-
∂
T
∂
n
.
Variation of local Nusselt number at hot wall and cold wall of cavity for different Rayleigh numbers is plotted in Figure 3. We calculate the average Nusselt number on the vertical boundary of the cavity at x
=
0 by
(83)
Nu
=
∫
0
1
‍
Nu
local
(
x
,
y
)
d
y
.
The average Nusselt numbers at the hot wall of the cavity are given in Table 7, which are according to the results of [16, 34, 35].
Figure4 presents streamlines and isotherms for Ra
=
1
0
4
,
1
0
5
,
1
0
6
,
1
0
7, respectively. For the streamline patterns, it is easy to see that elliptical vortex at the cavity center break up into two vortices tending to approach to the corners differentially heated sides of the cavity with Rayleigh number increasing. With the increase of Rayleigh numbers Ra, the temperature convection becomes increasingly prominent, the isotherms gradually transform into horizontal except for the immediate neighborhood of the hot and cold walls which remain parallel to the isothermal vertical walls. These results are keeping with the results of [16, 17, 34–36].Table 5
Comparison of maximum vertical velocity aty
=
0.5 with mesh size used in computation for Example 2.
Ra
Present
Reference [31]
Reference [34]
Reference [35]
Reference [16]
Reference [36]
1
0
4
19.6221 (64)
19.63 (64)
19.51 (41)
19.63 (71)
19.629 (32)
19.79 (101)
1
0
5
68.4791 (64)
68.48 (64)
68.22 (81)
68.85 (71)
68.65 (32)
70.63 (101)
1
0
6
220.444 (64)
220.46 (64)
216.75 (81)
221.6 (71)
220.57 (32)
227.11 (101)
1
0
7
695.416 (64)
694.14 (64)
—
702.3 (71)
699.3 (64)
714.48 (301)Table 6
Comparison of maximum horizontal velocity atx
=
0.5 with mesh size used in computation for Example 2.
Ra
Present
Reference [31]
Reference [34]
Reference [16]
Reference [36]
10
4
16.1825 (64)
16.19 (64)
16.18 (41)
16.183 (32)
16.10 (101)
10
5
34.7359 (64)
34.74 (64)
34.81 (81)
34.76 (32)
34 (101)
10
6
64.8122 (64)
64.81 (64)
65.33 (81)
64.81 (32)
65.40 (101)
10
7
146.914 (64)
148.40 (64)
—
148.8 (64)
143.56 (301)Table 7
Comparison of average Nusselt number on the vertical boundary of the cavity atx
=
0 with mesh size used in computation for Example 2.
Ra
Present
Reference [8]
Reference [34]
Reference [16]
Reference [36]
1
0
4
2.24511 (64)
2.21 (12)
2.24 (41)
2.245 (32)
2.254 (101)
1
0
5
4.52572 (64)
4.53 (22)
4.52 (81)
4.521 (32)
4.598 (101)
1
0
6
8.87446 (64)
9.00 (32)
8.92 (81)
8.826 (32)
8.976 (101)
1
0
7
16.9603 (64)
—
—
16.52 (64)
16.656 (301)Figure 1
Natural convection in cavity: the physical domain with its boundary conditions.Natural convection in cavity: comparison of vertical velocity at the midheight (a) and horizontal velocity at midwidth (b) for different Rayleigh numbers.
(a)
(b)Natural convection in cavity: comparison of local Nusselt numbers along the hot wall (x
=
0) (a) and the cold wall (x
=
1) (b) for different Rayleigh numbers.
(a)
(b)Natural convection in cavity: streamlines and isotherms forRa
=
1
0
4, 1
0
5, 1
0
6, 1
0
7 from left to right, respectively.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)


## 6. Conclusions

In this paper, a fully discrete stabilized finite element method based on the Crank-Nicolson linear extrapolation scheme in time for natural convection problem is given. The method is unconditionally stable. Optimal error estimates in H1 seminorm and suboptimal error estimates inL
2 norm are derived for velocity and temperature with a constrain condition. The derived theoretical results are supported by two numerical examples.



---

*Source: 393494-2014-05-25.xml*