A Finite Element Variational Multiscale Method Based on Two Local Gauss Integrations for Stationary Conduction-Convection Problems

Abstract

A new finite element variational multiscale (VMS) method based on two local Gauss integrations is proposed and analyzed for the stationary conduction-convection problems. The valuable feature of our method is that the action of stabilization operators can be performed locally at the element level with minimal additional cost. The theory analysis shows that our method is stable and has a good precision. Finally, the numerical test agrees completely with the theoretical expectations and the “ exact solution,” which show that our method is highly efficient for the stationary conduction-convection problems.

1. Introduction

The conduction-convection problems constitute an important system of equations in atmospheric dynamics and dissipative nonlinear system of equations. Many authors have worked on these problems [1–8]. The governing equations couple viscous incompressible flow and heat transfer process [9], where the incompressible fluid is the Boussinesq approximation to the nonstationary Navier-Stokes equations. Christon et al. [10] summarized some relevant results for the fluid dynamics of thermally driven cavity. A multigrid (MG) technique was applied for the conduction-convection problems [11, 12]. Luo et al. [13] combined proper orthogonal decomposition (POD) with the Petrov-Galerkin least squares mixed finite element (PLSMFE) method for the problems. In [14], a Newton iterative mixed finite element method for the stationary conduction-convection problems was shown by Si et al. In [15], Si and He gave a defect-correction mixed finite element method for the stationary conduction-convection problems. In [3], an analysis of conduction natural convection conjugate heat transfer in the gap between concentric cylinders under solar irradiation was carried out. In [16], Boland and Layton gave an error analysis for finite element methods for steady natural convection problems. Variational multiscale (VMS) method which defines the large scales in a different way, namely, by a projection into appropriate subspaces, see Guermond [17], Hughes et al. [18–20] and Layton [21], and other literatures on VMS methods [22–24]. The new finite element VMS strategy requires edge-based data structure and a subdivision of grids into patches. It does not require a specification of mesh-dependent parameters and edge-based data structure, and it is completely local at the element level. Consequently, the new VMS method under consideration can be integrated in existing codes with very little additional coding effort.

For the conduction-convection problems, we establish such system that Ω be a bounded domain in 𝑅^𝑑(𝑑=2or3), with Lipschitz-continuous boundary 𝜕Ω. In this paper, we consider the stationary conduction-convection problem as follows: −2𝜈∇⋅𝐷(𝑢)+(𝑢⋅∇𝑢)+∇𝑝=𝜆𝑗𝑇,𝑥∈Ω,∇⋅𝑢=0,𝑥∈Ω,−Δ𝑇+𝜆𝑢⋅∇𝑇=0,𝑥∈Ω,𝑢=0,𝑇=𝑇₀,𝑥∈𝜕Ω, (1.1) where 𝐷(𝑢)=(∇𝑢+∇𝑢^𝑇)/2 is the velocity deformation tensor, (𝑢,𝑝,𝑇)∈𝑋×𝑀×𝑊,Ω⊂𝑅^𝑑 is a bounded convex domain. 𝑢=(𝑢₁(𝑥),𝑢₂(𝑥))^𝑇 represents the velocity vector, 𝑝(𝑥) the pressure, 𝑇(𝑥) the temperature, 𝜆>0 the Grashoff number, 𝑗=(0,1)^𝑇 the two-dimensional vector and 𝜈>0 the viscosity.

The study is organized as follows. In the next section, the finite element VMS method is given. In Section 3, we give the stability. The error analysis is given in Section 4. In Section 5, we show some numerical test. The last but not least is the conclusion given in Section 6.

2. Finite Element VMS Method

Here, we introduce some notations 𝑋=𝐻¹₀(Ω)^𝑑,𝑀=𝐿²₀(Ω)=𝜑∈𝐿²(Ω);_Ω𝜑𝑑𝑥=0,𝑊=𝐻¹(Ω). (2.1) For ℎ>0, finite-dimension subspace (𝑋_ℎ,𝑀_ℎ,𝑊_ℎ)⊂(𝑋,𝑀,𝑊) is introduced which is associated with Ω_𝑒, a triangulation of Ω into triangles or quadrilaterals, assumed to be regular in the usual sense. In this study, the finite-element subspaces of personal preference are defined by setting the continuous piecewise (bi)linear velocity and pressure subspace, let 𝜏_ℎ be the regular triangulations or quadrilaterals of the domain Ω and define the mesh parameter ℎ=max_{Ω𝑒∈𝜏ℎ}{diam(Ω_𝑒)},𝑋_ℎ=𝑣∈𝑋∶𝑣|_Ω𝑒∈𝑅_𝑙Ω_𝑒^𝑑∀Ω_𝑒∈𝜏_ℎ,𝑀_ℎ=𝑞∈𝑀∶𝑞_Ω𝑒∈𝑅_𝑙Ω_𝑒∀Ω_𝑒∈𝜏_ℎ,𝑊_ℎ=𝜙∈𝑀∶𝜙_Ω𝑒∈𝑅_𝑙Ω_𝑒∀Ω_𝑒∈𝜏_ℎ, (2.2) where 𝑊_0ℎ=𝑊_ℎ∩𝐻¹₀, 𝑙≥1 is integers. 𝑅_𝑙(Ω_𝑒)=𝑃_𝑙(Ω_𝑒) if Ω_𝑒 is triangular and 𝑅_𝑙(Ω_𝑒)=𝑄_𝑙(Ω_𝑒) if Ω_𝑒 is quadrilateral. Here (𝑋_ℎ,𝑀_ℎ) does not satisfy the discrete Ladyzhenskaya-Babuška-Brezzi (LBB) condition sup_{𝑣ℎ∈𝑋ℎ}𝑑𝑣_ℎ,𝑝_ℎ∇𝑣_ℎ0≥𝛽𝑝_ℎ0,∀𝑝_ℎ∈𝑀_ℎ. (2.3) Now, in order to stabilize the convective term appropriately for the higher Reynolds number and avoid the extra storage, we supply finite element VMS method that the local stabilization form of the difference between a consistent and an underintegrated mass matrices based on two local Gauss integrations at element level as the stabilize term 𝐺𝑝_ℎ,𝑞_ℎ=𝜖_𝑑𝑎_𝑘𝑝_ℎ,𝑞_ℎ−𝑎₁𝑝_ℎ,𝑞_ℎ. (2.4) Here, 𝑎_𝑘𝑝_ℎ,𝑞_ℎ=𝑝^𝑇_𝐺𝑀_𝑘𝑞_𝐺,𝑎₁𝑝_ℎ,𝑞_ℎ=𝑝^𝑇_𝐺𝑀₁𝑞_𝐺,𝑝^𝑇_𝐺=𝑝₁,𝑝₂…,𝑝_𝑁^𝑇,𝑞_𝐺=𝑞₁,𝑞₂,…,𝑞_𝑁,𝑀_𝑖𝑗=𝜙_𝑖,𝜙_𝑗,𝑝_ℎ=^𝑁_𝑖=1𝑝_𝑖𝜙_𝑖,𝑝_𝑖=𝑝_ℎ𝑥_𝑖,∀𝑝_ℎ∈𝑀_ℎ𝑀,𝑖=1,2,…,𝑁,_𝑘=𝑀^𝑘_{𝑖𝑗𝑁×𝑁},𝑀₁=𝑀¹_{𝑖𝑗𝑁×𝑁}, (2.5) the stabilization parameter 𝜖_𝑑(𝜖_𝑑=𝑜(ℎ)) in this scheme acts only on the small scales, 𝜙_𝑖 is the basis function of the velocity on the domain Ω such that its value is one at node 𝑥_𝑖 and zero at other nodes, and 𝑁 is the dimension of 𝑀_ℎ. The symmetric and positive matrices 𝑀^𝑘_𝑖𝑗, 𝑘≥2 and 𝑀¹_𝑖𝑗 are the stiffness matrices computed by using 𝑘-order and 1-order Gauss integrations at element level, respectively. 𝑝_𝑖 and 𝑞_𝑖, 𝑖=1,2,…,𝑁 are the values of 𝑝_ℎ and 𝑞_ℎ at the node 𝑥_𝑖. In detail, the stabilized term can be rewritten as 𝐺𝑝_ℎ,𝑞_ℎ=𝜖_{𝑑Ω𝑒∈𝜏ℎΩ𝑒,𝑘}𝑝_ℎ𝑞_ℎ𝑑𝑥−_Ω𝑒,1𝑝_ℎ𝑞_ℎ𝑑𝑥,∀𝑝_ℎ,𝑞_ℎ∈𝑀_ℎ,𝐺𝑝,𝑞=𝑝−_ℎ𝑝,𝑞−_ℎ𝑞. (2.6)𝐿²-projection operator ∐_ℎ∶𝐿²(Ω)→𝑅₀ with the following properties [25]: 𝑝,𝑞_ℎ=∐_ℎ𝑝,𝑞_ℎ,∀𝑝∈𝑀,𝑞_ℎ∈𝑅₀;∐_ℎ𝑝₀≤𝑐𝑝₀∐,∀𝑝∈𝑀;𝑝−_ℎ𝑝₀≤𝑐ℎ𝑝₁,∀𝑝∈𝐻¹(Ω)∩𝑀. (2.7)

Lemma 2.1 (see [26]). Let (𝑋_ℎ,𝑀_ℎ) be defined as above, then there exists a positive constant 𝛽 independent of ℎ, such that 𝐵𝑢,𝑝;𝑣,𝑞≤𝑐‖𝑢‖₁+𝑝₀‖𝑣‖₁+𝑞₀,∈(),𝛽𝑢𝑢,𝑝𝑣,𝑞𝑋,𝑀_ℎ1+𝑣𝑝_ℎ0≤sup_{𝑣ℎ,𝑞ℎ∈𝑋ℎ,𝑀ℎ}𝐵𝑢_ℎ,𝑝_ℎ;𝑣_ℎ,𝑞_ℎ‖𝑣‖₁+𝑞₀,∀𝑢_ℎ,𝑝_ℎ∈𝑋_ℎ,𝑀_ℎ,𝐺𝑝,𝑞≤𝐶𝐼−𝐼𝐼_ℎ𝑝₀𝐼−𝐼𝐼_ℎ𝑞₀,∀𝑝,𝑞∈𝑀. (2.8)

Using the above notations, the VMS variational formulation of problems (1.1) reads as follows.

Find (A₁)(𝑢_ℎ,𝑝_ℎ,𝑇_ℎ)∈𝑋_ℎ×𝑀_ℎ×𝑊_ℎ such that 𝑎𝑢_ℎ,𝑣_ℎ−𝑑𝑝_ℎ,𝑣_ℎ+𝑑𝑞_ℎ,𝑢_ℎ+𝑏𝑢_ℎ,𝑢_ℎ,𝑣_ℎ+𝐺𝑝_ℎ,𝑞_ℎ=𝜆𝑗𝑇_ℎ,𝑣_ℎ,∀𝑣_ℎ∈𝑋_ℎ,𝜑_ℎ∈𝑀_ℎ;𝑎𝑇_ℎ,𝜓_ℎ+𝜆𝑏𝑢_ℎ,𝑇_ℎ,𝜓_ℎ=0,∀𝜓_ℎ∈𝑊_0ℎ. (2.9) Given (A₂)(𝑢_ℎ^𝑛−1,𝑇_ℎ^𝑛−1), find (𝑢^𝑛_ℎ,𝑝^𝑛_ℎ,𝑇^𝑛_ℎ)∈𝑋_ℎ×𝑀_ℎ×𝑊_ℎ such that 𝑎𝑢^𝑛_ℎ,𝑣_ℎ−𝑑𝑝^𝑛_ℎ,𝑣_ℎ+𝑑𝑞_ℎ,𝑢^𝑛_ℎ+𝑏𝑢^𝑛_ℎ,𝑢_ℎ^𝑛−1,𝑣_ℎ+𝑏𝑢_ℎ^𝑛−1,𝑢^𝑛_ℎ,𝑣_ℎ+𝐺𝑝^𝑛_ℎ,𝑞_ℎ=𝑏𝑢_ℎ^𝑛−1,𝑢_ℎ^𝑛−1,𝑣_ℎ+𝜆𝑗𝑇^𝑛_ℎ,𝑣_ℎ,∀𝑣_ℎ∈𝑋_ℎ,𝜑_ℎ∈𝑀_ℎ;𝑎𝑇^𝑛_ℎ,𝜓_ℎ+𝜆𝑏𝑢_ℎ^𝑛−1,𝑇^𝑛_ℎ,𝜓_ℎ=0,∀𝜓_ℎ∈𝑊_0ℎ, (2.10) where 𝑎(𝑢,𝑣)=𝜈(∇𝑢,∇𝑣),𝑎(𝑇,𝜓)=(∇𝑇,∇𝜓)𝑑(𝑞,𝑣)=(𝑞,div𝑣), and 𝑏()=))+1𝑢,𝑣,𝑤((𝑢⋅∇𝑣,𝑤2))=1((div𝑢𝑣,𝑤2))−1((𝑢⋅∇𝑣,𝑤2((𝑢⋅∇)𝑤,𝑣),𝑏=))+1𝑢,𝑇,𝜓((𝑢⋅∇𝑇,𝑤2()=1div𝑢𝑇,𝜓2()−1𝑢⋅∇𝑇,𝜓2(𝑢⋅∇)𝜓,𝑇. (2.11)(B₁) There exists a constant 𝐶 which only depends on Ω, such that(i)‖𝑢‖₀≤𝐶‖∇𝑢‖₀,‖𝑢‖_0,4≤𝐶‖∇𝑢‖₀,forall𝑢∈𝐻¹₀(Ω)^𝑑(or𝐻¹₀(Ω)),(ii)‖𝑢‖_0,4≤𝐶‖𝑢‖₁,forall𝑢∈𝐻¹(Ω)^𝑑(iii)‖𝑢‖_0,4≤2^1/2‖∇𝑢‖₀^1/2‖𝑢‖₀^1/2,forall𝑢∈𝐻¹₀(Ω)^𝑑(or𝐻¹₀(Ω))(B₂) Assuming 𝜕Ω∈𝐶^𝑘,𝛼(𝑘≥0,𝛼>0), then, for 𝑇₀∈𝐶^𝑘,𝛼(𝜕Ω), there exists an extension 𝑇₀ in 𝐶₀^𝑘,𝛼(𝑅^𝑑), such that 𝑇_0𝑘,𝑞≤𝜀,𝑘≥0,1≤𝑞≤∞, (2.12)where 𝜀 is an arbitrary positive constant.(B₃)𝑏(⋅,⋅,⋅)and𝑏(⋅,⋅,⋅) have the following properties.(i) For all 𝑢∈𝑋,𝑣,𝑤∈𝑋,𝑇⋅𝜑∈𝐻¹₀(Ω), there holds that 𝑏(𝑢,𝑣,𝑤)=−𝑏(𝑢,𝑤,𝑣),𝑏𝑢,𝑇,𝜓=−𝑏𝑢,𝜓,𝑇. (2.13)(ii) For all 𝑢∈𝑋,𝑣∈𝐻¹(Ω)^𝑑,𝑇∈𝐻¹(Ω),forall𝑤∈𝑋(or𝜑∈𝐻¹₀(Ω)), there holds that |𝑏(𝑢,𝑣,𝑤)|≤𝑁‖∇𝑢‖₀‖∇𝑣‖₀‖∇𝑤‖₀,𝑏𝑢,𝑇,𝜑≤𝑁‖∇𝑢‖₀‖∇𝑇‖₀∇𝜑₀, (2.14)where  𝑁=sup_{𝑢,𝑣,𝑤}|𝑏(𝑢,𝑣,𝑤)|/‖∇𝑢‖₀‖∇𝑣‖₀‖∇𝑤‖₀,𝑁=sup_{𝑢,𝑣,𝑤}|𝑏(𝑢,𝑇,𝜑)|/‖∇𝑢‖₀‖∇𝑇‖₀‖∇𝜑‖₀.

3. Stability Analysis

Lemma 3.1. The trilinear form 𝑏 satisfies the following estimate: 𝑏𝑢_ℎ,𝑣_ℎ,𝑤+𝑏𝑣_ℎ,𝑢_ℎ,𝑤+𝑏𝑤,𝑢_ℎ,𝑣_ℎ≤𝐶logℎ^1/2∇𝑢_ℎ0∇𝑣_ℎ0‖𝑤‖₀. (3.1)

Theorem 3.2. Suppose that (B₁)-(B₃) are valid and 𝜀 is a positive constant number, such that 64𝐶²𝑁𝜆𝜀3𝜈²<1,16𝐶²𝜆²𝑁𝜀3𝜈<1,∇𝑇₀₀≤𝜀4,𝑇₀₀≤𝐶𝜀4. (3.2) Then (𝑢^𝑚_ℎ,𝑇^𝑚_ℎ) defined by (A₂) satisfies ∇𝑢^𝑚_ℎ0≤8𝐶²𝜆𝜀,3𝜈∇𝑇^𝑚_ℎ0≤𝜀. (3.3)

Proof. We prove this theorem by the inductive method. For 𝑚=1, (3.3) holds obviously. Assuming that (3.3) holds for 𝑚=𝑛−1, we want to prove that it holds for 𝑚=𝑛. We estimate ‖Δ𝑢^𝑛_ℎ‖ firstly. Letting 𝑣_ℎ=𝑢^𝑛_ℎ,𝑞_ℎ=0 in the first equation of (2.10) and using (2.13), we get 𝑎𝑢^𝑛_ℎ,𝑢^𝑛_ℎ+𝑏𝑢^𝑛_ℎ,𝑢_ℎ^𝑛−1,𝑢^𝑛_ℎ=𝑏𝑢_ℎ^𝑛−1,𝑢_ℎ^𝑛−1,𝑢^𝑛_ℎ+𝜆𝑗𝑇^𝑛_ℎ,𝑢^𝑛_ℎ. (3.4) Setting 𝑇_ℎ^𝑛−1=𝑘_ℎ^𝑛−1+𝑇₀ and using (2.14), we have 𝜈∇𝑢^𝑛_ℎ0≤𝑁∇𝑢^𝑛_ℎ0∇𝑢_ℎ^𝑛−1₀+𝑁∇𝑢_ℎ^𝑛−12₀+𝐶²𝜆∇𝑘_ℎ^𝑛−1₀+𝐶𝜆∇𝑇₀₀. (3.5) Letting 𝑇^𝑛_ℎ=𝑘^𝑛_ℎ+𝑇₀,𝜓=𝑘^𝑛_ℎ in the second equation of (2.10), we can obtain 𝑎𝑘^𝑛_ℎ,𝑘^𝑛_ℎ=−𝜆𝑏𝑢^𝑛_ℎ,𝑇₀,𝑘^𝑛_ℎ−𝑎𝑇₀,𝑘^𝑛_ℎ. (3.6) Using (2.12), (2.14), and (3.2), we get ∇𝑘_ℎ^𝑛−1₀≤𝜆𝑁∇𝑢_ℎ^𝑛−1₀∇𝑇₀₀+∇𝑇₀₀≤𝜆𝑁𝜀4∇𝑢_ℎ^𝑛−1₀+∇𝑇₀₀≤3𝜀8≤3𝜀4,𝜈−𝑁∇𝑢_ℎ^𝑛−1₀∇𝑢^𝑛_ℎ0≤𝑁∇𝑢_ℎ^𝑛−12₀+𝐶²𝜆𝜀≤𝐶²𝜆𝜀+64𝐶⁴𝑁9𝜈²𝜆²𝜀²≤4𝐶²𝜆𝜀3. (3.7) Using (3.2), we have 𝜈−𝑁‖∇𝑢_ℎ^𝑛−1‖₀≥7𝜈/8. Then, ∇𝑢^𝑛_ℎ0≤8𝐶²𝜆𝜀.3𝜈 (3.8) Combining (2.12), (2.14), (3.2), and (3.6), we arrive at ∇𝑘^𝑛_ℎ0≤𝜆𝑁∇𝑢^𝑛_ℎ0∇𝑇₀₀+∇𝑇₀₀≤3𝜀4,∇𝑇^𝑛_ℎ0≤∇𝑘^𝑛_ℎ0+∇𝑇₀₀≤𝜀. (3.9)(3.10) Therefore, we finish the proof.

4. Error Analysis

In this section, we establish the 𝐻¹-bound of the error 𝑢^𝑛_ℎ−𝑢,𝑇^𝑛_ℎ−𝑇 and 𝐿²-bounds of the error 𝑝^𝑛_ℎ−𝑝. Setting (𝑒^𝑛,𝜇^𝑛,𝜂^𝑛)=(𝑢^𝑛_ℎ−𝑢_ℎ,𝑝^𝑛_ℎ−𝑝_ℎ,𝑇^𝑛_ℎ−𝑇_ℎ). Firstly, we give some Lemmas.

Lemma 4.1. In [4], If B₁-B₃ hold, (𝑢,𝑝, 𝑇)∈𝐻^𝑚+1(Ω)×𝐻^𝑚(Ω)×𝐻^𝑚+1(Ω) and (𝑢_ℎ,𝑝_ℎ,𝑇_ℎ)∈𝑋_ℎ×𝑀_ℎ×𝑊_ℎ are the solution of problem (A₁) and (A₂), respectively, then there holds that ∇(𝑢−𝑢_ℎ)₀+𝑝−𝑝_ℎ0+∇(𝑇−𝑇_ℎ)₀≤Ch^𝑚‖𝑢‖_𝑚+1+𝑝_𝑚+‖𝑇‖_𝑚+1. (4.1)

Lemma 4.2. Under the assumptions of Theorem 3.2, (A₂) has a unique solution (𝑢_ℎ,𝑝_ℎ,𝑇_ℎ)∈ 𝑋_ℎ×𝑀_ℎ×𝑊_ℎ, such that 𝑇|_𝜕Ω=𝑇₀ and ∇𝑢_ℎ0≤8𝐶²𝜆𝜀,3𝜈∇𝑇_ℎ0≤𝜀. (4.2)

The detail proof we can see [4, 13, 14].

Theorem 4.3. Under the assumption of Theorem 3.2, there holds ∇𝑒^𝑛₀≤𝐶²𝜆𝜀2^𝑛−3,3𝜈∇𝜂^𝑛₀≤𝜀2^𝑛+1,𝜇^𝑛₀≤𝛽⁻¹𝜈𝜀2+4𝐶²𝜆𝜀3()𝐶,𝑛=1𝜈+2𝑁𝜀²𝜆𝜀2^𝑛−3𝐶3𝜈+𝑁²𝜆𝜀2^𝑛−43𝜈²+𝐶²𝜆𝜀2^𝑛,𝑛≥2. (4.3)

Proof. Subtracting (2.10) from (2.9), we get the following error equations, namely (𝑒^𝑛,𝜇^𝑛,𝜂^𝑛) satisfies 𝑎𝑒^𝑛,𝑣_ℎ−𝑑𝜇^𝑛,𝑣_ℎ+𝑑𝑞_ℎ,𝑒^𝑛+𝑏𝑒^𝑛,𝑢_ℎ^𝑛−1,𝑣_ℎ+𝑏𝑢_ℎ^𝑛−1,𝑒^𝑛,𝑣_ℎ+𝐺𝜇^𝑛,𝑞_ℎ=𝑏𝑒^𝑛−1,𝑒^𝑛−1,𝑣_ℎ+𝜆𝑗𝜂^𝑛,𝑣_ℎ,𝑎𝜂^𝑛,𝜓_ℎ+𝜆𝑏𝑒^𝑛,𝑇^𝑛_ℎ,𝜓_ℎ+𝜆𝑏𝑢_ℎ^𝑛−1,𝜂^𝑛,𝑞_ℎ=0. (4.4)(4.5) Here, let 𝜓_ℎ=𝜂^𝑛, in (4.5), then we have 𝑎𝜂^𝑛,𝜂^𝑛+𝜆𝑏𝑒^𝑛,𝑇^𝑛_ℎ,𝜂^𝑛=0. (4.6) By using (2.14), we get ∇𝜂^𝑛₀≤𝜆𝑁𝜀∇𝑒^𝑛₀. (4.7) In (4.4), we take 𝑣_ℎ=𝑒^𝑛∈𝑋_ℎ,𝑞_ℎ=𝜇^𝑛, then 𝑎𝑒^𝑛,𝑒^𝑛+𝑏𝑒^𝑛,𝑢_ℎ^𝑛−1,𝑒^𝑛+𝑏𝑢_ℎ^𝑛−1,𝑒^𝑛,𝑒^𝑛+𝐺𝜇^𝑛,𝜇^𝑛=𝑏𝑒^𝑛−1,𝑒^𝑛−1,𝑒^𝑛+𝜆𝑗𝜂^𝑛,𝑒^𝑛. (4.8) Using (2.13) and (2.14), we have 𝜈∇𝑒^𝑛₀+𝐺𝜇^𝑛,𝜇^𝑛≤𝑁∇𝑒^𝑛₀∇𝑢_ℎ^𝑛−1₀+𝑁∇𝑒^𝑛−12₀+𝐶²𝜆∇𝜂^𝑛₀, (4.9) then, we obtain 𝜈−𝑁∇𝑢_ℎ^𝑛−1₀∇𝑒^𝑛₀≤𝑁∇𝑒^𝑛−12₀+𝐶²𝜆∇𝜂^𝑛₀. (4.10) By using 𝜈−𝑁‖∇𝑢_ℎ^𝑛−1‖₀≥7𝜈/8. Equations (3.3) and (4.2), we get 78𝜈∇𝑒^𝑛₀≤𝑁∇𝑒^𝑛−12₀+𝐶²𝜆∇𝜂^𝑛₀≤𝑁∇𝑢_ℎ^𝑛−1₀+𝑁∇𝑢_ℎ0+𝐶²𝜆²𝑁𝜀∇𝑒^𝑛−1₀≤16𝑁𝐶²𝜆𝜀3𝜈+𝐶²𝜆²𝑁𝜀∇𝑒^𝑛−1₀=7𝜈16∇𝑒^𝑛−1₀,∇𝑒^𝑛₀≤12∇𝑒^𝑛−1₀. (4.11) From the inductive method, we know, for 𝑛=1, subtracting (2.10) from (2.9), we can get 𝑎𝑒¹,𝑣_ℎ−𝑑𝜇¹,𝑣_ℎ+𝑑𝑞_ℎ,𝑒¹+𝑏𝑢_ℎ,𝑢_ℎ,𝑣_ℎ+𝐺𝜇¹,𝑞_ℎ=𝜆𝑗𝑇^𝑛,𝑣_ℎ. (4.12) Letting 𝑣_ℎ=𝑒¹,𝑞_ℎ=𝜇¹ in (4.12) and using (2.14), we have ∇𝑒¹₀+𝐺𝜇¹,𝜇¹≤𝜈⁻¹𝑁∇𝑢_ℎ²₀+𝜈⁻¹𝐶²𝜆∇𝑇_ℎ0≤64𝐶⁴𝜆²𝑁𝜀²9𝜈³+𝐶²𝜆𝜀𝜈≤4𝐶²𝜆𝜀,3𝜈 (4.13) then ∇𝑒¹₀≤4𝐶²𝜆𝜀.3𝜈 (4.14) By (4.7), we have ∇𝜂¹₀≤𝜆𝑁𝜀∇𝑒^𝑛₀≤𝜆𝑁4𝐶²𝜆𝜀²≤𝜀3𝜈4. (4.15) Letting 𝑞_ℎ=0 in (4.12), (2.14), and (3.9), using Lemma 2.1, we get 𝛽𝜇¹₀≤𝜈∇𝑒¹₀+𝑁∇𝑢_ℎ²₀+𝐶𝜆𝑇_ℎ0≤𝜈𝜀2+4𝐶²𝜆𝜀.3𝜈 (4.16) Assuming that (4.3) is true for 𝑛=𝑘−1, using (4.7) and (4.11), we know that both of them are valid for 𝑛=𝑘. Using (4.7) holds for 𝑛=𝑘, we let 𝑞_ℎ=0 in (4.4) and using Lemma 2.1, (4.5), and (3.3), we have 𝛽𝜇^𝑛₀≤(𝜈+2𝑁𝜀)∇𝑒^𝑛₀+𝑁∇𝑒^𝑛−12₀+𝐶²𝜆∇𝜂^𝑛−1₀≤()𝐶𝜈+2𝑁𝜀²𝜆𝜀2^𝑛−3𝐶3𝜈+𝑁²𝜆𝜀2^𝑛−43𝜈²+𝐶²𝜆𝜀2^𝑛. (4.17)

Theorem 4.4. Under the assumptions of Theorem 4.3, then there holds that lim_𝑛→∞𝑢^𝑛_ℎ−𝑢_ℎ^𝑛−1₀+∇𝑢^𝑛_ℎ−𝑢_ℎ^𝑛−1₀=0,∇𝑒^𝑛₀+𝜇^𝑛₀+∇𝜂^𝑛₀≤𝐹logℎ^1/2∇𝑢^𝑛_ℎ−𝑢_ℎ^𝑛−1₀𝑢^𝑛_ℎ−𝑢_ℎ^𝑛−1₀+𝐻𝜀2^𝑛+1, (4.18) where 𝐹 and 𝐻 are two positive constants.

Proof. By using (B₁) and triangle inequality, we have 𝑢^𝑛_ℎ−𝑢_ℎ^𝑛−1₀+∇(𝑢^𝑛_ℎ−𝑢_ℎ^𝑛−1)₀≤(𝐶+1)∇𝑒^𝑛₀+∇𝑒^𝑛−1₀. (4.19) Using Theorem 4.3, letting 𝑛→∞, we obtain (4.18). Taking 𝑣_ℎ=𝑒^𝑛,𝑞_ℎ=𝜇^𝑛 in (4.4) and using (2.14), we get 𝑎𝑒^𝑛,𝑒^𝑛+𝑏𝑒^𝑛,𝑢_ℎ^𝑛−1,𝑒^𝑛+𝐺𝜇^𝑛,𝜇^𝑛=−𝑏𝑢^𝑛_ℎ−𝑢_ℎ^𝑛−1,𝑢^𝑛_ℎ−𝑢_ℎ^𝑛−1,𝑒^𝑛+𝜆𝑗𝜂^𝑛,𝑒^𝑛. (4.20) By (2.14) and Lemma 3.1, we deduce 𝜈−𝑁𝑢_ℎ^𝑛−1₀∇𝑒^𝑛₀+𝐺𝜇^𝑛,𝜇^𝑛≤𝐹logℎ^1/2∇𝑢^𝑛_ℎ−𝑢_ℎ^𝑛−1₀𝑢^𝑛_ℎ−𝑢_ℎ^𝑛−1₀+𝐹²𝜆∇𝜂^𝑛₀. (4.21) Combining (3.3) and (4.7), we obtain 𝜈−8𝑁𝜀3𝜈∇𝑒^𝑛₀≤𝐹logℎ^1/2∇𝑢^𝑛_ℎ−𝑢_ℎ^𝑛−1₀𝑢^𝑛_ℎ−𝑢_ℎ^𝑛−1₀+𝐹²𝜆²𝑁𝜀∇𝑒^𝑛−1₀. (4.22) Using (3.2), we get ∇𝑒^𝑛₀≤𝐹logℎ^1/2∇(𝑢^𝑛_ℎ−𝑢_ℎ^𝑛−1)₀𝑢^𝑛_ℎ−𝑢_ℎ^𝑛−1₀+𝐻𝜀2^𝑛+1. (4.23) Combining (3.2), (4.7), and (4.17), we get ∇𝜂^𝑛₀≤𝐹logℎ^1/2∇(𝑢^𝑛_ℎ−𝑢_ℎ^𝑛−1)₀𝑢^𝑛_ℎ−𝑢_ℎ^𝑛−1₀+𝐻𝜀2^𝑛+1,∇𝜇^𝑛₀≤𝐹logℎ^1/2∇(𝑢^𝑛_ℎ−𝑢_ℎ^𝑛−1)₀𝑢^𝑛_ℎ−𝑢_ℎ^𝑛−1₀+𝐻𝜀2^𝑛+1. (4.24) Here, we complete the proof.

Theorem 4.5. Under the assumptions of Theorem 4.3, the following inequality: ∇(𝑢−𝑢^𝑛_ℎ)₀+𝑝−𝑝^𝑛_ℎ0+∇(𝑇−𝑇^𝑛_ℎ)₀≤𝐹₁ℎ^𝑚‖𝑢‖_𝑚+1+𝑝_𝑚+‖𝑇‖_𝑚+1+𝐹logℎ^1/2∇𝑢^𝑛_ℎ−𝑢_ℎ^𝑛−1₀𝑢^𝑛_ℎ−𝑢_ℎ^𝑛−1₀+𝐻𝜀2^𝑛+1, (4.25) holds, where 𝐹₁ and 𝐻 are the positive constants.

Proof. By Lemma 4.1, Theorem 4.4, and the triangle inequality, this theorem is obviously true.

5. Numerical Test

This section presents the numerical results that complement the theoretical analysis.

5.1. Convergence Analysis

In our experiment, Ω=[0,1]×[0,1] is the unit square in 𝑅². Let 𝑇₀=0 on left and lower boundary of the cavity, 𝜕𝑇/𝜕𝑛=0 on upper boundary of the cavity, and 𝑇₀=4𝑦(1−𝑦) on right boundary of the cavity (see Figure 1). Physics model of the cavity flows: 𝑡=0, that is, 𝑛=0 initial values on boundary. In general, we cannot know the exact solution of the stationary conduction-convection equations. In order to get the exact solution, we design the procedure as follows. Firstly, solving the stationary conduction-convection equations by using the 𝑃₂-𝑃₁-𝑃₂ finite element pair, which holds stability, on the finer mesh, we take the solution as the exact solution. Secondly, the absolute error is obtained by comparing the exact solution and the finite element solutions with VMS methods. Finally, we can easily obtain errors and convergence rates.

(a)

(b)

(c)

Figure 1 

From (a) to (c): physics model of the cavity flows, vertical midlines for Re = 2000, ℎ = 1/100, horizontal midlines for Re = 2000, ℎ = 1/100.

5.2. Driven Cavity

In this experiment, Ω=[0,1]×[0,1] is the unit square in 𝑅². Let 𝑇₀=0 on left and lower boundary of the cavity, 𝜕𝑇/𝜕𝑛=0 on upper boundary of the cavity, and 𝑇₀=4𝑦(1−𝑦) on right boundary of the cavity (see Figure 1). Physics model of the cavity flows: 𝑡=0, that is, 𝑛=0 initial values on boundary. Solving the stationary conduction-convection equations by using the 𝑃₂-𝑃₁-𝑃₂ finite element pair, which holds stability results, on the finer mesh, we take the solution as the exact solution. From Figures 1 and 2, we know that the solution of finite element VMS using 𝑃₁-𝑃₁-𝑃₁ element agree completely with the “exact solution.” In Figure 3, we choose Re=2000, divide the cavity into 𝑀×𝑁=100×100, from left to right shows the numerical streamline, the numerical isobar, and the numerical isotherms. In Figure 4, we choose Re=3000, divide the cavity into 𝑀×𝑁=100×100, from left to right shows the numerical streamline, the numerical isobar, and the numerical isotherms.

(a)

(b)

Figure 2 

From (a) to (b): vertical midlines for Re = 3000, ℎ = 1/100, horizontal midlines for Re = 3000, ℎ = 1/100.

(a)

(b)

(c)

Figure 3 

For Re = 2000, ℎ = 1/100, from (a) to (c): velocity streamlines, the pressure level lines, numerical isotherms.

(a)

(b)

(c)

Figure 4 

For Re = 3000, ℎ = 1/100, from (a) to (c): velocity streamlines, the pressure level lines, numerical isotherms.

Remark 5.1. Our VMS finite element method based on two local Gauss integrations and 𝜀_𝑑=0.1ℎ is suitable for the Sobolev space. Throughout the paper, our analysis and numerical tests are all carried out for the 𝑃₁-𝑃₁-𝑃₁ element (see Tables 1 and 2).

Table 1 

VMS: 𝑃₁-𝑃₁-𝑃₁ element.

Table 2 

VMS: 𝑃₁-𝑃₁-𝑃₁ element.

6. Conclusion

In this paper, we studied a finite element VMS algorithm based on two local Gauss integrations to solve the stationary conduction-convection problem. From Figures 1 and 2, we see that the solution of VMS using 𝑃₁-𝑃₁-𝑃₁ and 𝜀_𝑑=0.1ℎ agrees completely with the “exact solution,” which shows that our method is highly efficient for the stationary conduction-convection problems. Numerical tests tell us that VMS finite element method based on two local Gauss integrations is very effective.

Acknowledgments

The project is supported by NSF of China (10971164) and the Research Foundation of Xianyang Normal University (06xsyk265).