An Adaptive Finite Element Method for Stationary Incompressible Thermal Flow Based on Projection Error Estimation

Abstract

An adaptive finite element method is presented for the stationary incompressible thermal flow problems. A reliable a posteriori error estimator based on a projection operator is proposed and it can be computed easily and implemented in parallel. Finally, three numerical examples are given to illustrate the efficiency of the adaptive finite element method. We also show that the adaptive strategy is effective to detect local singularities in the physical model of square cavity stationary flow in the third example.

1. Introduction

As we all know the stationary incompressible thermal flow occurs in many fields of the industrial sectors and natural world. The mathematical model of this flow can be described by a set of coupled equations, which consist of the conservation of mass, momentum, and energy equations. The stationary incompressible thermal flow problems constitute an important system of equations in atmospheric dynamics and dissipative nonlinear system of equations. The main difficulties for the numerical simulation of the stationary incompressible thermal flow problem include not only the incompressibility and strong nonlinearity, but also the the coupling between the energy equation and the equations governing the fluid motion. And the numerical methods for solving these problems have attracted the attention of many researchers. Many authors have worked hard to study for a great variety of efficient numerical schemes for these equations [1–6]. In [2], Si et al. gave a coupled Newton iterative mixed finite element method for stationary conduction convection problems. In [4], a nonconforming mixed finite element method for the stationary conduction convection problems was shown by Shi and Ren. In [6], a new finite element variational multiscale (VMS) method based on two local Gauss integrations was proposed and analyzed for the stationary conduction convection problems by Jiang et al. Furthermore, there are also some scholars devoted to the nonstationary incompressible thermal flow problems [7–11]. In [10], a reduced mixed finite element formulation based on proper orthogonal decomposition for the nonstationary conduction convection problems is given by Luo and his coworkers. In [11], Li et al. gave a fully discrete FVM formulation for the nonstationary conduction convection problems and derived the error analysis between the fully discrete FVM solution and the accurate solution.

In the numerical simulation of stationary incompressible thermal flow problems, still a big challenge is how to increase the accuracy of the numerical approximations for the solutions. The overall accuracy of numerical approximations often deteriorates due to local singularities or large variations in small scales. Motivated by this phenomenon, a direct strategy is that the grids near the critical regions are refined adaptively to improve the quality of the approximate solutions. In 1978, Babuska and Rheinboldt developed a mathematical theory for a class of a posteriori error estimator of finite element solutions [12, 13]. Based on a posteriori error estimator, the adaptive finite element methods have been designed to generate optimal or near optimal meshes and compute solutions more exactly to the problems with boundary layers. Many researchers are dedicated to deriving such a posteriori error estimators and have obtained many good results in recent decades [14–16]. There are a lot of work devoted to the development of the a posteriori analysis [17–20]. In 1996, Verfürth built a basic foundation for a posteriori error estimation for nonlinear equations [18], based on which Zhang et al. [21] derived a residual a posteriori error estimator for finite element approximations of stationary incompressible thermal flow problems, in 2011. In this paper, a posteriori error estimator is presented based on a projection operator. Its actions can use only standard nodal data structures, so it can be computed locally at the element level and implemented in parallel. The global upper bound for the error of the finite element discretization is derived following some given assumptions.

This paper is organized as follows. In Section 2, we introduce the governing equations, the notations and some preliminary results used for the stationary incompressible thermal flow problems. The a posteriori error estimator based on local projection operator is presented in Section 3. In Section 4, three numerical results and the numerical analysis to validate the effectiveness of the adaptive method are laid out. The first two examples are known solution examples, and the last example is a physical model of square cavity stationary flow. Finally, we give a short conclusion in Section 5.

2. Governing Equations and the Functional Setting

In this section, we consider the stationary incompressible thermal flow problems in two-dimensions. The viscous incompressible flow and heat transfer for the fluid satisfy the incompressible Navier-Stokes equations coupled with the energy conservation equation under the Boussinesq hypothesis is as follows.

Find 𝐮=(𝑢₁,𝑢₂),𝑝,𝑇 such that ()()−𝜈Δ𝐮+𝐮⋅∇𝐮+∇𝑝=𝜆𝑗𝑇,∀𝑥∈Ω,div𝐮=0,∀𝑥∈Ω,−Δ𝑇+𝜆𝐮⋅∇𝑇=0,∀𝑥∈Ω,𝐮=0,𝑇=𝑇₀,∀𝑥∈Γ, (1) where Ω is a bounded domain in 𝑅² assumed to have Lipschitz continuous boundary Γ. 𝐮=(𝑢₁,𝑢₂) represents the velocity vector, 𝑝 the pressure, 𝑇 the temperature, 𝜆>0 the Groshoff number, 𝑗=(0,1) the unit vector, and 𝜈>0 the viscosity.

We introduce the following Hilbert spaces 𝐗=𝐻¹₀(Ω)², 𝐕={𝐯∈𝐗,∇⋅𝐯=0}, 𝑌=𝐿²₀(Ω)={𝑞∈𝐿²(Ω),∫_Ω𝑞d𝑥=0}, and 𝑊=𝐻¹(Ω), 𝑊₀=𝐻¹₀(Ω). The standard variational formulation of (1) is given by: find (𝐮,𝑝,𝑇)∈𝐗×𝑌×𝑊 and 𝑇|_Γ=𝑇₀ satisfying 𝜈𝑎₁(𝐮,𝐯)−𝑑𝑝,𝐯+𝑏₁()𝑑𝑎𝐮,𝐮,𝐯=𝜆𝑗𝑇,𝑣∀𝐯∈𝐗,𝑞,𝐮=0∀𝑞∈𝑌,₂𝑇,𝜓+𝜆𝑏₂𝐮,𝑇,𝜓=0∀𝜓∈𝑊₀,𝐮=0,𝑇=𝑇₀∀𝑥∈Γ. (2) Here we used the notations 𝑎₁()=()𝑎𝐮,𝐯∇𝐮,∇𝐯,∀𝐮,𝐯∈𝐗,₂𝑇,𝜓=∇𝑇,∇𝜓,∀𝑇∈𝑊,𝜓∈𝑊₀,𝑑=𝑏𝑞,𝐯𝑞,div𝐯,∀𝐯∈𝐗,𝑞∈𝑌,₁()=1𝐮,𝐯,𝐰2_Ω²_𝑖,𝑘=1𝑢_𝑖𝜕𝑣_𝑘𝜕𝑥_𝑖𝑤_𝑘−d𝑥_2𝑖,𝑘=1𝑢_𝑖𝜕𝑤_𝑘𝜕𝑥_𝑖𝑣_𝑘𝑏d𝑥,∀𝐮,𝐯,𝐰∈𝐗,₂=1𝐮,𝑇,𝜓2_Ω²_𝑖=1𝑢_𝑖𝜕𝑇𝜕𝑥_𝑖𝜓d𝑥−_Ω²_𝑖=1𝑢_𝑖𝜕𝜓𝜕𝑥_𝑖𝑇d𝑥,∀𝐮∈𝐗,𝑇∈𝑊,𝜓∈𝑊₀, (3) and (⋅,⋅) the inner product in 𝐿²(Ω) or in its vector value versions. The norm and seminorm in 𝐻^𝑘(Ω)² are denoted by ‖⋅‖_𝑘 and |⋅|_𝑘, respectively. 𝐻¹₀(Ω) will denote the closure of 𝐶^∞₀(Ω) with respect to the norm ‖⋅‖₁. The space 𝐗 is equipped with the norm ‖∇⋅‖₀ or its equivalent norm ‖⋅‖₁ due to the Poincaré inequality. Throughout this paper, we use the letter 𝐶 (with or without subscripts) to denote a generic positive constant which may stand for different value at its different occurrences.

For the finite element discretization, let 𝐾_ℎ be the regular triangulations of the domain Ω, indexed by a parameter ℎ=max_{𝐾∈𝐾ℎ}{ℎ_𝐾;ℎ_𝐾=diam(𝐾)}. We choose the finite element subspace 𝐗_𝐡⊆𝐗,𝑌_ℎ⊆𝑌,𝑊_ℎ⊆𝑊 as follows: 𝐗_𝐡=𝐯_𝐡∈𝑐(Ω)²;𝐯_𝐡|_𝐾∈𝑃₂(𝐾)²,∀𝐾∈𝐾_ℎ,𝑌_ℎ=𝑞_ℎ∈𝑐(Ω);𝑞_ℎ|_𝐾∈𝑃₁(𝐾),∀𝐾∈𝐾_ℎ,𝑊_ℎ=𝜓_ℎ∈𝑐(Ω);𝜓_ℎ|_𝐾∈𝑃₂(𝐾),∀𝐾∈𝐾_ℎ, (4) where 𝑃_𝑙(𝐾)(𝑙=1,2) is the space of piecewise polynomials of degree 𝑙 on 𝐾. We will also need the piecewise constant space 𝑅₀=𝑙_ℎ∈𝐿²(Ω),𝑙_ℎ|_𝐾∈𝑃₀(𝐾),∀𝐾∈𝐾_ℎ. (5) We denote that 𝑊_0ℎ=𝑊_ℎ⋂𝐻¹₀(Ω). It is obvious that the (𝐗_𝐡,𝑌_ℎ) satisfies the discrete LBB condition sup_{𝐯𝐡∈𝐗𝐡}𝑑𝐯_𝐡,𝑞_ℎ∇𝐯_𝐡0≥𝛽𝑞_ℎ0∀𝑞_ℎ∈𝑌_ℎ. (6)

With the previous notations, the Galerkin finite element discretization of (2) is given by the following: find (𝐮_𝐡,𝑝_ℎ,𝑇_ℎ)∈𝐗_𝐡×𝑌_ℎ×𝑊_ℎ such that 𝜈𝑎₁𝐮_𝐡,𝐯_𝐡−𝑑𝑝_ℎ,𝐯_𝐡+𝑏₁𝐮_𝐡,𝐮_𝐡,𝐯_𝐡=𝜆𝑗𝑇_ℎ,𝐯_𝐡,∀𝐯_𝐡∈𝐗_𝐡,𝑑𝑞_ℎ,𝐮_𝐡=0,∀𝑞_ℎ∈𝑌_ℎ,𝑎₂𝑇_ℎ,𝜓_ℎ+𝜆𝑏₂𝐮_𝐡,𝑇_ℎ,𝜓_ℎ=0,∀𝜓_ℎ∈𝑊_0ℎ. (7)

For problem (1) the following assumptions and results are recalled (see [8, 22, 23]).(𝐀_𝟏) There exists a constant 𝐶₀ which only depends on Ω such that (i)‖𝐮‖₀≤𝐶₀‖∇𝐮‖₀,‖𝐮‖_0,4≤𝐶₀‖∇𝐮‖₀,∀𝐮∈𝐻¹₀(Ω)²,(ii)‖𝐮‖_0,4≤𝐶₀‖𝐮‖₁,∀𝐮∈𝐻¹(Ω)²,(iii)‖𝐮‖_0,4≤2^1/2‖∇𝐮‖₀^1/2‖𝐮‖₀^1/2,∀𝐮∈𝐻¹₀(Ω)².(𝐀_𝟐) Assuming 𝜕Ω∈𝐶^𝑘,𝛼(𝑘≥0,𝛼≥0), then for 𝑇₀∈𝐶^𝑘,𝛼(𝜕Ω) and a sufficient small 𝜀, there exists an extension in 𝐶₀^𝑘,𝛼(𝑅²) (denoted 𝑇₀ also), such that ‖𝑇₀‖_𝑘,𝑞≤𝜀,𝑘≥0,1≤𝑞≤∞, where 𝜀 is an arbitrary small positive constant number.(𝐀_𝟑)𝑏₁(⋅,⋅,⋅) and 𝑏₂(⋅,⋅,⋅) have the following properties: (i) for all𝐮,𝐯, 𝐰∈𝐗 (or 𝑇,𝜓∈𝐻¹₀(Ω)), there holds that 𝑏₁(𝐮,𝐯,𝐯)=0,𝑏₁(𝐮,𝐯,𝐰)=−𝑏₁(),𝑏𝐮,𝐰,𝐯₂(𝐮,𝑇,𝑇)=0,𝑏₂𝐮,𝑇,𝜓=−𝑏₂𝐮,𝜓,𝑇, (8)(ii) for all 𝐮∈𝐗, 𝐯∈𝐻¹(Ω)², 𝐰∈𝐗, (or T ∈𝐻¹(Ω), and 𝜓∈𝐻¹₀(Ω)), there holds that 𝑏₁(𝐮,𝐯,𝐰)≤𝑁₁‖∇𝐮‖₀‖∇𝐯‖₀‖∇𝐰‖₀,𝑏₂𝐮,𝑇,𝜓≤𝑁₂‖∇𝐮‖₀‖∇𝑇‖₀∇𝜓₀, (9)where𝑁₁=sup_{𝐮,𝐯,𝐰}𝑏₁(𝐮,𝐯,𝐰)‖∇𝐮‖₀‖∇𝐯‖₀‖∇𝐰‖₀,𝑁₂=sup_{𝐮,𝑇,𝜓}𝑏₂𝐮,𝑇,𝜓‖∇𝐮‖₀‖∇𝑇‖₀∇𝜓₀. (10)

Theorem 1 (see [8]). Under the assumption of (𝐀_𝟏)~(𝐀_𝟑), and letting 𝐴=2𝜈⁻¹𝜆(3𝐶₀+1)‖𝑇₀‖₁, 𝐵=2‖∇𝑇₀‖₀+2(𝐶₀𝜆)⁻¹𝐴, assuming there exist 0<𝜎₁,𝜎₂≤1 such that 𝜈⁻¹𝑁𝐴≤1−𝜎₁, 𝜎₁⁻¹𝜈⁻¹𝐶²₀𝜆²𝐵𝑁≤1−𝜎₂, then there exists a unique weak solution (𝐮,𝑝,𝑇)∈𝐗×𝑌×𝑊, for problem (1) and ‖∇𝐮‖₀≤𝐴,‖∇𝑇‖₀≤𝐵. (11)

Theorem 2 (see [8]). Under the assumption of (𝐀_𝟏)~(𝐀_𝟑) and 64𝐶²₀𝑁𝜆𝜀/3𝜈²<1, 16𝐶²₀𝜆²𝑁𝜀/3𝜈<1, ‖∇𝑇₀‖₀≤𝜀/4, ‖𝑇₀‖₀≤𝐶₀𝜀/4, then (7) has a unique solution (𝐮_𝐡,𝑝_ℎ,𝑇_ℎ)∈𝐗_𝐡×𝑌_ℎ×𝑊_ℎ such that 𝑇|_𝜕Ω=𝑇₀ and ∇𝐮_𝐡0≤8𝐶²₀𝜆3𝜈𝜀,∇𝑇_ℎ0≤𝜀. (12)

3. A Projection Error Estimation for Stationary Incompressible Thermal Flow Problems

Before deriving the projection estimator, we define the orthogonal projection operator [24] Π∶𝐿²(Ω)→𝑅₀ acting on the pressure which satisfies the following properties: (𝐼−Π)𝑔,𝑙_ℎ=0,∀𝑔∈𝐿²(Ω),𝑙_ℎ∈𝑅₀,Π𝑔₀≤𝐶𝑔₀,∀𝑔∈𝐿²(Ω),(𝐼−Π)𝑔₀≤𝐶ℎ𝑔₁,∀𝑔∈𝐻¹(Ω). (13) Here, 𝐼 is the identity operator. In the following discussion, operator Π^2×2, Π^2×1, which acts on the velocity deformation tensor and the gradient of the temperature ∇𝑇_ℎ, are also denoted by Π for simplicity.

Now, based on the residual between the gradient of the finite element solutions velocity component ∇𝐮_𝐡; pressure component 𝑝_ℎ; temperature component ∇𝑇_ℎ; and their projections Π∇𝐮_𝐡, Π∇𝑝_ℎ, and Π∇𝑇_ℎ; our projection estimator can be constructed locally as follows: 𝜂_Π,𝐾∶=(𝐼−Π)∇𝐮_𝐡0,𝐾+(𝐼−Π)𝑝_ℎ0,𝐾+(𝐼−Π)∇𝑇_ℎ0,𝐾. (14)

Then the global error estimator is given by 𝜂∶=_{𝐾∈𝐾ℎ}𝜂²_Π,𝐾^1/2. (15)

Remark 3. Based on orthogonal projection properties of operator Π, the computation of 𝐼−Π can be done at the element level using only standard nodal data structures.

Remark 4. Our local projection error estimator can be computed more precisely and explicitly based on two local Gauss integrations technique presented in [25, 26]. We give the detailed form of the 𝐼−Π operator as folows: (𝐼−Π)𝑞_ℎ0,𝐾=_{𝐾,𝑗(𝑗≥2)}𝑞_ℎ²d𝑥−_𝐾,1𝑞_ℎ²d𝑥^1/2,∀𝑞_ℎ∈𝑌_ℎ. (16) Here∫_𝐾,𝑙𝑔(𝑥)d𝑥,𝑙=1,𝑗 denotes a quadrature formula.
Before giving the global upper bound, we recall a lemma in [27].

Lemma 5. There exists a positive constant C such that 𝐶ℎ∇𝑝_ℎ0≤(𝐼−Π)𝑝_ℎ0∀𝑝_ℎ∈𝑌_ℎ. (17)

Remark 6. This lemma was given in [27], the proof of this lemma will be rewritten here.

Proof. We note that 𝑝_ℎ is continuous and Π𝑝_ℎ is a constant on each element 𝐾 from the definition of space 𝑌_ℎ and the projection operator Π, so it is obvious that ∇(Π𝑝_ℎ)|_𝐾=0.
Additionally, under some mild assumption on 𝐾_ℎ, the following inverse inequalities holds [28]. ∇𝑞_ℎ0≤𝐶_𝐼ℎ⁻¹𝑞_ℎ0,𝑞_ℎ∈𝑌_ℎ. (18) Spaces concerning with vector-valued functions will hold too.
As a result, using the previous inverse inequality (18), ℎ²∇𝑝_ℎ²₀=_{𝐾∈𝐾ℎ}ℎ²∇𝑝_ℎ²_0,𝐾=_{𝐾∈𝐾ℎ}ℎ²∇(𝐼−Π)𝑝_ℎ²_0,𝐾≤_{𝐾∈𝐾ℎ}𝐶²_𝐼(𝐼−Π)𝑝_ℎ²_0,𝐾≤𝐶(𝐼−Π)𝑝_ℎ²₀. (19)
Then, this completes the proof of Lemma 5.

For ∇𝐮_𝐡∈𝑅₁^2×2, ∇𝑇_ℎ∈𝑅₁^2×1, and spaces consisting of vector-valued functions will also hold for 𝐶ℎ∇𝐮_𝐡1≤(𝐼−Π)∇𝐮_𝐡0,∀𝐮_𝐡∈𝐗_𝐡,𝐶ℎ∇𝑇_ℎ1≤(𝐼−Π)∇𝑇_ℎ0,∀𝑇_ℎ∈𝑊_ℎ. (20)

Theorem 7. There is a constant C depending on the smallest angle in the triangulation 𝐾_ℎ and the domain Ω, such that, we have the following global upper bounds: ∇(𝐮−𝐮_𝐡)₀+𝑝−𝑝_ℎ0+∇(𝑇−𝑇_ℎ)₀≤𝐶𝜂, (21) where C is independent of h.

Proof. We assume that there are two positive constants 𝐶_𝑚,𝐶_𝑛, such that 𝐶_𝑚≤𝑝₁,𝑝_ℎ1,‖𝐮‖₂,∇𝐮_𝐡1,‖𝑇‖₂,∇𝑇_ℎ1≤𝐶_𝑛 (22) we can deduce the following results with the discussion of Lemma 5: ∇(𝐮−𝐮_𝐡)₀+𝑝−𝑝_ℎ0+∇(𝑇−𝑇_ℎ)₀≤𝐶ℎ‖𝐮‖₂+𝐶ℎ𝑝₁+𝐶ℎ‖𝑇‖₂𝐶≤𝐶^𝑛𝐶_𝑚ℎ∇𝐮_𝐡1𝐶+𝐶^𝑛𝐶_𝑚ℎ𝑝_ℎ1𝐶+𝐶^𝑛𝐶_𝑚ℎ∇𝑇_ℎ1≤𝐶(𝐼−Π)∇𝐮_𝐡0+(𝐼−Π)𝑝_ℎ0+(𝐼−Π)∇𝑇_ℎ0≤𝐶𝜂. (23)

Then, the theorem is proved.

4. Numerical Experiments

In this section, we present three examples to illustrate the effectiveness of the adaptive method. In all experiments, the nonlinear systems are solved by the Newton iteration [1] and implemented in the two-dimensional framework using a public domain finite element software [29]. The first two examples are known solution problems. The last example deals with “lid driven cavity,” which is a benchmark problem for testing numerical schemes.

We make use of the Newton iteration scheme presented in [1] as follows. Given (𝐮_𝐡^𝐧−𝟏,𝑇_ℎ^𝑛−1), find (𝐮^𝐧_𝐡,𝑝^𝑛_ℎ,𝑇^𝑛_ℎ)∈𝑋_ℎ×𝑌_ℎ×𝑊_ℎ and 𝑇^𝑛_ℎ|_Γ=𝑇₀ such that 𝜈𝑎₁𝐮^𝐧_𝐡,𝐯_𝐡−𝑑𝑝^𝑛_ℎ,𝐯_𝐡+𝑑𝑞_ℎ,𝐮^𝐧_𝐡+𝑏₁𝐮^𝐧_𝐡,𝐮_𝐡^𝐧−𝟏,𝐯_𝐡+𝑏₁𝐮_𝐡^𝐧−𝟏,𝐮^𝐧_𝐡,𝐯_𝐡=𝑏₁𝐮_𝐡^𝐧−𝟏,𝐮_𝐡^𝐧−𝟏,𝐯_𝐡+𝜆𝑗𝑇_ℎ^𝑛−1,𝐯_𝐡,∀𝐯_𝐡∈𝐗_𝐡,∀𝑞_ℎ∈𝑌_ℎ,𝑎₂𝑇^𝑛_ℎ,𝜓_ℎ+𝜆𝑏₂𝐮^𝐧_𝐡,𝑇^𝑛_ℎ,𝜓_ℎ=0,∀𝜓_ℎ∈𝑊_0ℎ. (24)

For the sake of clarity, the main idea is offered about the refinement strategy presented in [29] for the algorithm of the discrete problem (7). From the beginning of original triangulation 𝐾⁰_ℎ, a polygonal approximation Ω and a series of refined triangulations 𝐾^𝑗_ℎ were constructed as below. Given 𝐾^𝑗_ℎ, first of all, we compute the error estimator 𝜂^𝑗_Π,𝐾 after computing the solution (𝐮_𝐣,𝑝_𝑗,𝑇_𝑗) from the problem (7); then the new mesh size is given by the following formulation: ℎ_𝑗+1=ℎ_𝑗𝑓_𝑗𝜂^𝑗_Π,𝐾, (25) here ℎ_𝑗 is the previous “meshsize” field, and 𝑓_𝑗 is a user function defined by 𝑓_𝑗𝜂=minmax^𝑗Π,𝐾𝑐𝜂^𝑗,1.0000,3, (26) where 𝜂^𝑗 means the mean of 𝜂^𝑗_Π,𝐾 and 𝑐 is an user coefficient generally close to one and the numbers 1.0000,3 also can be changed by the user according to the requirement.

For simplicity, we choose 𝑐=0.8 in all later numerical tests. Certainly, we can also consider other refinement strategy, the one in [17] for example, whereas there is no much differences between them.

The computation adaptive strategy here is to choose a tolerance 𝜂^∗, start from the original triangulation 𝐾⁰_ℎ, and then compute the 𝜂 by the following two steps.

Step 1. If 𝜂≤𝜂^∗, stop the calculation, then we obtain the final solution; otherwise, go to the Step 2.

Step 2. Compute the 𝜂_Π,𝐾 and their mean value 𝜂, and generate a new mesh size ℎ by the above strategy and then solve the problem (7) and recompute 𝜂 based on this new triangulation. Go back to Step 1.

For convenience of presentation, we introduce the following notation: (i)𝑁∶= number of elements in 𝐾^𝑗_ℎ,(ii)Eff∶=𝜂/err the effective index, that is, the ratio between the error estimator and the true error. Here, 𝜂²∶=∑_{𝐾∈𝐾ℎ}𝜂²_Π,𝐾, err∶={𝑒₁(𝐮)²+𝑒₀(𝑝)²+𝑒₁(𝑇)²}^1/2,(iii)𝑒₁(𝐮)∶=‖∇(𝐮−𝐮_𝐡)‖₀, 𝑒₀(𝑝)∶=‖𝑝−𝑝_ℎ‖₀, 𝑒₁(𝑇)∶=‖∇(𝑇−𝑇_ℎ)‖₀.

Example 1. For this problem the analytical solutions were taken to be 𝑢₁𝑥,𝑦=10𝑥²(𝑥−1)²𝑦,𝑢𝑦−12𝑦−1₂𝑥,𝑦=−10𝑥(𝑥−1)(2𝑥−1)𝑦²𝑦−1²,(),(𝑎)𝑝=102𝑥−12𝑦−1𝑇=cosh−cosh𝑎𝑦[(𝑎)],cosh−1 (27) with the chosen function added to the right-hand side of (1). The problem is solved on the unite square. From Figure 1, we can easily find that the temperature field has large variation near the boundary 𝑦=1, when 𝑎=50. So we expect grid points are clustered in the region of rapid variation of the temperature field to improve accuracy. In terms of the meshes refined on the left of Figure 2, we notice that the adaptive strategy creates a lot of triangles in the area near the boundary of 𝑦=1. We also give the numerical solution for temperature field on the right of Figure 2. It shows very good agreement with our predictions.
Tables 1–4 present the numerical solutions based on uniform meshes and the adaptive procedure. As shown in Tables 1 and 2, we can obviously observe that the adaptive strategy almost gets the same good results compared with the uniform one when 𝑎=10 in which case the temperature field shows smoothly and flat. But uniform refinements do not achieve good approximation when setting 𝑎=50, where the temperature field has a severe variation in small scales, yet the adaptive strategy still obtains much better approximation solution. For example, when the global total error decreases to 0.10297, the uniform meshes needed increase to 8712, while only 1049 triangles are used to get a smaller global total error in the adaptive procedure by comparison of Tables 3 and 4.
All in all, our projection a posteriori error estimator can assess the domain of rapid solution variation and refine them to improve the global accuracy with less meshes.

Table 1 

Uniform meshes for 𝜈=1.0; 𝜆=1.0; 𝑎=10.

Table 2 

A series of adapted meshes for 𝜈=1.0; 𝜆=1.0; 𝑎=10.

Table 3 

Uniform meshes for 𝜈=1.0,𝜆=1.0,𝑎=50.

Table 4 

A series of adapted meshes for 𝜈=1.0,𝜆=1.0,𝑎=50.

(a)

(b)

Figure 1 

Exact solution for temperature. 𝑎=10 (a) and 𝑎=50 (b).

(a)

(b)

Figure 2 

Final adapted mesh (a) and the numerical solution for temperature (b) for 𝑎=50.

Example 2. We choose the known solutions as follows, Ω=[0,1]×[0,1], and the chosen functions are added to the right-hand side of (1) such that the exact solution of the problem is given by 𝑇𝑥,𝑦=𝑢₁𝑥,𝑦+𝑢₂,𝐮=𝑢𝑥,𝑦𝑥,𝑦₁𝑥,𝑦,𝑢₂,𝑢𝑥,𝑦₁=𝑒𝑥,𝑦1−cos2𝜋^𝑟1𝑥−1𝑒^𝑟1𝑒−1×sin2𝜋^𝑟2𝑦−1𝑒^𝑟2𝑟−1²𝑒2𝜋(𝑒^{𝑟2𝑦𝑟2}),𝑢−1₂𝑒𝑥,𝑦=−sin2𝜋^𝑟1𝑥−1𝑒^𝑟1×𝑒−11−cos2𝜋^𝑟2𝑦−1𝑒^𝑟2𝑟−1¹𝑒2𝜋(𝑒^{𝑟1𝑥𝑟1}),𝑝−1𝑥,𝑦=𝑟₁𝑟₂𝑒sin2𝜋^𝑟1𝑥−1𝑒^𝑟1𝑒−1×sin2𝜋^𝑟2𝑦−1𝑒^𝑟2𝑒−1^𝑟1𝑥𝑒^𝑟2𝑦(𝑒^𝑟1−1)(𝑒^𝑟2),−1 (28) where 𝑟₁ and 𝑟₂ are two strictly positive real parameters. The velocity field of this solution is similar to a counter clockwise vortex in a unite box. we can move the center of this vortex that has coordinates 𝑥₀=(1/𝑟₁)log((𝑒^𝑟1+1)/2) and 𝑦₀=(1/𝑟₂)log((𝑒^𝑟2+1)/2) along with the change of the parameters 𝑟₁ and 𝑟₂. The center goes rapidly towards the right-hand vertical side when increasing 𝑟₁, in the same way it approaches the top edge when increasing 𝑟₂. We can verify clearly this phenomenon from the right of Figures 3 and 4. We can also observe that the pressure has large variation in small scale near the right-hand vertical side and the top edge from the left of Figures 3 and 4. So we expect that the adaptive strategy can be sensitive to detect these critical regions and refine them to avoid the numerical oscillation.
Figure 5 presents a sequence of adaptive meshes for 𝑟₁=8.5 and 𝑟₂=0.1. The mesh refinement appears near the right-hand vertical side after a series of adaptive iterations. Also the adaptive mesh refinement for 𝑟₁=3.5 and 𝑟₂=9.1 in Figure 6 appears near the corner of the right-top, where the solutions have high gradients. These results are consistent with our analysis of problem's requirement. Finally, we give the numerical solutions after a series of adapted iterations when 𝑟₁=8.5,𝑟₂=0.1, and 𝑟₁=3.5,𝑟₂=9.1 in Figures 7 and 8.
To show the effectiveness of our adaptive method and the adaptive procedures based on the residual posteriori error estimator [21, 30, 31], in which the detailed formulation of the local error estimator is presented.
We present the numerical results for Example 2 in Tables 5–10. Tables 5 and 6 show the results on uniform meshes and adaptive meshes when 𝑟₁=8.5 and 𝑟₂=0.1. The results on adaptive meshes based on residual a posteriori error estimator are presented in Table 7. Also, we report the same numerical results in Tables 8–10, when 𝑟₁=3.5 and 𝑟₂=9.1. From Tables 5–10 we notice that the effective index Eff confirms the reliability and efficiency of posteriori error indicator 𝜂.
Comparing Tables 5 and 6, we observe that the global error of the adaptive procedures decreases much faster than those obtained by the quasiuniform ones. For example, when the error is around 0.18, we need 11858 triangles by uniform procedures, while only 2074 triangles are needed by adaptive ones when 𝑟₁=8.5 and 𝑟₂=0.1. For another example, when the error decreases to 0.13428 with 2168 meshes by adaptive procedures; however, the global error by uniform ones only drops to 0.173739 even if the meshes required increases to 12800 from the comparison of Tables 8 and 9, when 𝑟₁=3.5and𝑟₂=9.1. Concerning the meshes, we find that our method cost less memory space than common Galerkin finite element method. In fact, it can nearly save four fifths the memory space of the computer to get almost the same precision.

Table 5 

Uniform meshes for 𝜈=1.0,𝜆=1.0,𝑟₁=8.5,𝑟₂=0.1.

Table 6 

A successive sequence of adapted meshes for 𝜈=1.0,𝜆=1.0,𝑟₁=8.5,𝑟₂=0.1.

Table 7 

A successive sequence of adapted meshes based on a residual error estimation for 𝜈=1.0,𝜆=1.0,𝑟₁=8.5,𝑟₂=0.1.

Table 8 

Uniform meshes for 𝜈=1.0,𝜆=1.0,𝑟₁=3.5, and 𝑟₂=9.1.

Table 9 

A successive sequence of adaptive meshes for 𝜈=1.0,𝜆=1.0,𝑟₁=3.5, and 𝑟₂=9.1.

Table 10 

A successive sequence of adaptive meshes based on a residual error estimation for 𝜈=1.0,𝜆=1.0, and 𝑟₁=3.5, 𝑟₂=9.1.

(a)

(b)

Figure 3 

Exact solution for 𝑟1=8.5and𝑟2=0.1; pressure (a) and velocity field (b).

(a)

(b)

Figure 4 

Exact solution for 𝑟1=3.5and𝑟2=9.1; pressure (a) and velocity field (b).

(a)

(b)

(c)

(d)

(e)

(f)

Figure 5 

The successive adaptive refined meshes for 𝑟1=8.5 and 𝑟2=0.1.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 6 

The successive adaptive refined meshes for 𝑟1=3.5 and 𝑟2=9.1.

(a)

(b)

Figure 7 

Numerical solution for 𝑟1=8.5 and 𝑟2=0.1; pressure (a) and velocity field (b).

(a)

(b)

Figure 8 

Numerical solution for 𝑟1=3.5 and 𝑟2=9.1; pressure (a) and velocity field (b).

Comparing Tables 6 and 7, we notice that we need 6994 triangles by our adaptive method less than 7468 triangles based on posteriori error estimator when the global error is about 0.062 with 𝑟₁=8.5,𝑟₂=0.1. We also notice that our adaptive iterations needed are less such that we can save more CPU time to obtain nearly the same accuracy. This means that our adaptive method are more sensitive to those critical regions and refine them more efficiently than the adaptive procedures based on the residual a posteriori error estimator. The numerical tests in Tables 9 and 10 also verify this conclusion with the global error of about 0.06 when 𝑟₁=3.5 and 𝑟₂=9.1.

Combining all discussions above, we can derive a conclusion that these results suffice to show that our adaptive strategy obtains much better approximation with less meshes and CPU time. That is to say we can save lots of work by our adaptive procedures.

Example 3 (Square cavity stationary thermal flow). The last example is a physical model of square cavity stationary flow, which is a popular benchmark problem for testing numerical schemes. The side length of the square cavity and the boundary conditions are given in Figure 9. From Figure 9, we can see that 𝑇=0 on left and lower boundaries, 𝜕𝑇/𝜕𝑛=0 on upper boundary, and 𝑇=4𝑦(1−𝑦) on right boundary of the cavity. In this example, we set 𝜆=1.0 and 𝜈=1.0.
The left part of Figure 10 is the initial mesh. Then the successive adaptive refined meshes are generated on the basis of a posteriori error estimator (14). The first and third adaptive meshes are presented in the right of Figure 10 and the left of the Figure 11, respectively. From these adaptively generated meshes, we observe that the adaptive strategy is capable of recognizing the singularities and the regions with high gradients of the solutions. In addition, we show the numerical solution of (𝐮_𝐡,𝑝_ℎ,𝑇_ℎ) after three levels of adaptive meshes refinement. The streamline of velocity numerical solutions appears on the right of Figure 11. The distributions of the pressure and temperature are provided on the left and right of Figure 12, respectively.

Figure 9 

Physical model of the cavity flows.

(a)

(b)

Figure 10 

Initial mesh (a). The first adaptation mesh (b).

(a)

(b)

Figure 11 

The third adaptation mesh (a). The streamline of velocity numerical solutions after three levels of adaptation mesh refinement (b).

(a)

(b)

Figure 12 

Numerical isobar of pressure solution (a) and numerical isotherm of temperature solutions (b) after three levels of adaptation mesh refinement.

In general, we cannot know the exact solutions of the stationary incompressible thermal flow problems. We take the numerical solutions using the 𝑃2−𝑃1−𝑃2 finite element pairs on the finer mesh 150×150 as the “exact solutions.” In Figure 13, we give the numerical solutions on the final adaptive meshes 3631 and uniform meshes 45000 concerning the vertical midline for the 𝑥 component of velocity and the horizontal midline for 𝑦 component of the velocity. We can see that the adaptive finite element method approach the standard Galerkin method on the much fine mesh.

(a)

(b)

Figure 13 

Vertical midline for 𝑥 component of velocity (a) and horizontal midline for 𝑦 component of velocity (b).

5. Conclusion

In this paper, we present an adaptive finite element method based on a posteriori error estimator for the stationary incompressible thermal flow problems. This error estimator is constructed by a projection operator and can be calculated explicitly and precisely by the difference of two Gauss integrations technique. The discussions and the numerical tests indicate that the projection error estimator is effective and efficient. Whereas there are still many questions in further analysis, such as the convergence and optimality of the adaptive finite element methods.

Acknowledgments

This work is supported by the National Natural Science Foundation of China (nos. 11001216, 11371288, 11371289, 11272251, and 11271298) and the China Scholarship Council (no. 201206285018).