A projection-based stabilized finite element method for steady-state natural convection problem

J. Math. Anal. Appl. 381 (2011) 469–484

Contents lists available at ScienceDirect

Journal of Mathematical Analysis and Applications www.elsevier.com/locate/jmaa

A projection-based stabilized finite element method for steady-state natural convection problem Aytekin Çıbık a , Songül Kaya b,∗ a b

Department of Mathematics, Gazi University, 06500, Ankara, Turkey Department of Mathematics, Middle East Technical University and Institute of Applied Mathematics, 06531, Ankara, Turkey

a r t i c l e

i n f o

a b s t r a c t

Article history: Received 16 December 2010 Available online 16 February 2011 Submitted by Goong Chen

We formulate a projection-based stabilization finite element technique for solving steadystate natural convection problems. In particular, we consider heat transport through combined solid and fluid media. This stabilization does not act on the large flow structures. Based on the projection stabilization idea, finite element error analysis of the problem is investigated and optimal errors for the velocity, temperature and pressure are established. We also present some numerical tests which both verify the theoretical predictions and demonstrate the method’s promise. © 2011 Elsevier Inc. All rights reserved.

Keywords: Projection-based method Finite element method Error analysis Natural convection equation

1. Introduction In this article we consider a projection-based numerical stabilization for a convection dominated coupled problem. This type of stabilization, which adds an eddy viscosity stabilization only on the fine scales, is introduced in a variationally consistent way for the stationary convection diffusion problem in [18] and for the Navier–Stokes problem in [13]. This report studies an extension of the projection-based subgrid stabilization finite element method for the steady-state natural convection problem. The steady-state natural convection problem including solid media in dimensionless form is given by

− Pr u + (u · ∇)u + ∇ p = Pr Ra T e in Ω f , ∇ · u = 0 in Ω f , u = 0 on ∂Ω f ,

u≡0

in Ω − Ω f = Ωs ,

−∇ · (κ ∇ T ) + (u · ∇) T = γ in Ω, ∂T T = 0 on Γ T , = 0 on Γ B , ∂n

(1.1)

for the velocity u, the pressure p and the temperature T in a regular bounded open domain Ω ⊂ Rd (d = 2, 3), with disjoint polyhedral domains Ωs , Ω f . Here Γ T = ∂Ω\Γ B where Γ B is a regular open subset of ∂Ω , γ is a forcing function, e is a unit vector in the direction of gravitational acceleration and Pr, Ra, κ > 0 refer to the Prandtl, Rayleigh numbers and thermal conductivity parameter, respectively. Furthermore, we consider the case κ ≡ κ f in Ω f and κ ≡ κs in Ωs where κ f and κs are positive constants. The system (1.1) uses Boussinesq approximation as governing equations.

*

Corresponding author. E-mail addresses: [email protected] (A. Çıbık), [email protected] (S. Kaya).

0022-247X/$ – see front matter doi:10.1016/j.jmaa.2011.02.020

©

2011 Elsevier Inc. All rights reserved.

470

A. Çıbık, S. Kaya / J. Math. Anal. Appl. 381 (2011) 469–484

The natural convection equation (1.1), which is also known as buoyancy-driven flows, above occurs in many practical problems. Some of the commonly used buoyancy driven flows are observed in nature; such as atmospheric fronts, katabatic winds etc., and in industry; such as dense gas dispersion, natural ventilation, solar collectors, insulation with double pane window, cooling of electronic equipments, cooling of nuclear reactors, etc.. Classical natural convection problem in fluid mechanics occurs in an enclosed domain [16]. For natural convection in enclosures, a boundary layer forms near the walls. Outside this layer, a rolling core is formed inside the enclosure. The boundary layer and the core could not be considered independent since the core is covered by the layer. There is a coupling between the core and the boundary layer. This coupling is the main reason of the difficulty in solving these systems analytically. Thus, numerical methods and experimental analysis are used. System (1.1) presents severe computational problems for large Rayleigh numbers. It is well known that, the solution of (1.1) is unique under some restrictions on the Rayleigh and Prandtl numbers. Uniqueness is lost for high Rayleigh numbers [24]. Standard Galerkin finite element method for natural convection problem usually yields inaccurate approximate solutions and may exhibit global spurious oscillations [7,23]. This disappointing behavior occurs since such methods lose stability and cannot adequately approximate solutions inside layers due to the dominance of convection terms and the strong coupling between the velocity, pressure and temperature. Various stabilization techniques of finite element methods with more satisfactory performance have been developed in [25]. Our work is directed towards the efficient stabilization technique of natural convection problem to avoid some drawbacks of the classical methods. Overviews of some common stabilization mechanisms for convection diffusion equation and Oseen problem were given in [4,22]. One type of the stabilization mechanisms is the projection-based stabilization [12,8,18,13]. The philosophy of the projection-based stabilization is to use projections into appropriate function spaces in order to decompose solution scales. In this way, the stabilization is added in different ways. A noteworthy Guermond’s stabilization idea of subgrid viscosity concept makes the diffusion acts only on the finest resolved mesh scale [8], with the definition of solution spaces via bubble functions. Based on the ideas developed in [12,8], several multiscale decompositions have been proposed in the literature [18,13,14]. Since then, considerable progress has been made for the use of projection-based stabilization method both in mathematical and computational analysis in past years [15,11]. The objective of this paper is to provide finite element error analysis of the projection-based stabilization method for solving steady-state natural convection equations. In the meantime, some numerical analysis and numerical results for the time-dependent natural convection equation can be found in literature [19,2]. To do authors best knowledge, the error estimates of the projection-based stabilization is applied to steady-state natural convection equations are not yet available. In this paper, we consider the same type projection-based stabilization technique of the steady-state Navier Stokes equations [15]. As in [15], we also define the large scale spaces on a coarser grid for the solution scales. Main difference in the present work comes from the technical point of view, which is the coupling of the Navier–Stokes equation to the energy equation. We first present stabilized finite element scheme and give comprehensive error analysis of this coupled problem. We derive error estimations for the velocity, temperature and pressure and show that these errors are optimal with respect to the mesh sizes along with the choices of viscosity parameters. To evaluate the performance and accuracy of the method, we provide numerical experiments. The plan of this paper is as follows. In Section 2, the variational formulation of the problem is derived and the projectionbased finite element scheme is presented for the steady-state natural convection problem. Notations and mathematical preliminaries are given in Section 3. Existence, uniqueness and stability properties of the discrete problems are given in Section 4. Section 5 contains error estimations for the velocity and temperature. Section 6 is devoted to the error estimation of the pressure. Section 7 includes two numerical experiments: one is standard benchmark problem of buoyancy-driven flow in an enclosed domain. The next numerical experiment is chosen to illustrate the convergence theorem. Conclusion follows in Section 8. 2. Scheme The following well-known functional vector spaces are considered to define a variational formulation of (1.1).





X := H10 (Ω f ) = u ∈ H1 (Ω f ): u = 0 on ∂Ω f ,





W := S ∈ H 1 (Ω): S = 0 on Γ B ,





Q := p ∈ L 2 (Ω):



p dx = 0 , Ω

V := H10,div (Ω f ) = {u ∈ X: ∇ · u = 0 in Ω f },

(·,·) denotes the L 2 (Ω) inner product. We remark that the vector-valued functions are denoted with boldface character. We introduce the following bilinear and trilinear forms, for u, v, w ∈ X , T , S ∈ W and q ∈ Q :

A. Çıbık, S. Kaya / J. Math. Anal. Appl. 381 (2011) 469–484

471

 a0 (u, v) =

∇ u : ∇ v dx,

(2.1)

κ ∇ T · ∇ S dx,

(2.2)

Ωf



a1 ( T , S ) = Ω



b(v, q) = −

q∇ · v dx,

Ωf

c 0 (u, w, v) =



1



2

(2.3)

 (u · ∇)v · w − (u · ∇)w · v dx,

(2.4)

 (u · ∇) T S − (u · ∇) S T dx,

(2.5)

Ωf

c 1 (u, T , S ) =



1



2 Ωf

 d ( T , v) =

T e · v dx.

(2.6)

Ωf

The variational formulation of (1.1) reads as follows: seek u ∈ X , p ∈ Q , T ∈ W such that

Pr a0 (u, v) + c 0 (u, u, v) + b(v, p ) = Pr Ra d( T , v), b(u, q) = 0, a1 ( T , S ) + c 1 (u, T , S ) = (γ , S )

(2.7)

for all (v, q, S ) ∈ ( X , Q , W ). The notations in Eqs. (2.7) are inspired by the work in [3], in which the standard Galerkin finite element method for (2.7) is studied. The scheme introduces the addition of global stabilization and then subtracts its effect onto large scales of the coupled equations for both velocity and temperature spaces. In this way, stabilization acts only on the smallest resolved scales of both scales. Let F H , G K be a conforming triangulation of Ω and let F h , G k be a refinement of F H , G K , i.e. H  h and K  k respectively. Let X h ⊂ X , W k ⊂ W and Q h ⊂ Q be conforming finite element spaces satisfying the discrete inf–sup condition (3.2) in Section 3 and L H , M K denote the finite element subspaces of ( L 2 (Ω))d . The discretization we investigate adds additional diffusion acting on all discrete velocity and temperature scales and then anti-diffuses on the scales resolvable on F H , G K as follows: find uh ∈ X h , p h ∈ Q h , T k ∈ W k , F H ∈ L H and G K ∈ M K such that





Pr a0 uh , vh +



b uh , q



 h

= 0,

H















F − ∇ uh , l H = 0,





        a1 T , S + α2 ∇ T k − G K , ∇ S k + c 1 uh , T k , S k = γ , S k ,  K  G − ∇ T k , m K = 0, k

k







α1 ∇ uh − F H , ∇ vh + c0 uh , uh , vh + b vh , ph = Pr Ra d T k , vh ,



(2.8)

(2.9) (2.10) (2.11)

for all (vh , qh , l H , S k , m K ) ∈ ( X h , Q h , L H , W k , M K ) where α1 := α1 (h) and α2 := α2 (k) are non-negative constant functions and user selected stabilization parameters. These parameters can be thought of as an additional viscosity in the coarse space. Remark 2.1. Multiscale decomposition requires selection of large scale spaces for both velocity and temperature, L H and M K , respectively. If both of them are selected as zero subspaces, then Galerkin formulation is recovered in [3]. We employ L H = ∇ X H and M K = ∇ W K choices of [18] for the large scale spaces to obtain the bounds in this paper. Some other possible choices for these spaces are L H ⊆ ∇ X h and M K ⊆ ∇ W k (see [14]). Let V h = {vh ∈ X H : (qh , ∇ · vh ) = 0, for all qh ∈ Q h } be the space of discretely divergence free functions. It is easy to verify the following: (2.9) and (2.11) imply that F H and G K are L 2 projections of ∇ uh and ∇ T k onto L H and M K , respectively. If we denote these projections with P H and P K , respectively, the properties of the projection operator give the reformulations of (2.8)–(2.11) in V h as follows: find uh ∈ V h , T k ∈ W k such that

















A 0 uh , vh + c 0 uh , uh , vh = Pr Ra d T k , vh ,



k

A1 T , S

k

h

k

+ c1 u , T , S

k





= γ,S

k



(2.12) (2.13)

472

A. Çıbık, S. Kaya / J. Math. Anal. Appl. 381 (2011) 469–484

for all (vh , S k ) ∈ ( V h , W k ) where













A 0 uh , vh = Pr a0 uh , vh + α1 ( I − P H )∇ uh , ( I − P H )∇ vh ,



k

A1 T , S

k





k

= a1 T , S

k





k

k

(2.14)



+ α2 ( I − P K )∇ T , ( I − P K )∇ S .

(2.15)

3. Notation and preliminaries We present some function spaces and their norms. The standard notations in [1] for Sobolev and Lebesgue spaces are used. The inner product and norm in ( L 2 (Ω))d , d = 2, 3 are denoted by (·,·) and  · . The norm in ( H k (Ω))d is denoted by  · k and the norms in Lebesgue spaces ( L p (Ω))d , 1  p < ∞, p = 2 by  ·  L p . For vanishing boundary values, we define H 01 (Ω) and its dual space, H −1 (Ω) and its norm is defined by

f−1 = sup v∈ X

|(f, v)| ∇ v

where (·,·) denotes the duality pairing. We make use of well-known Sobolev embedding theorem for the following spaces: if Ω is bounded and has a Lipschitz boundary then H 1 (Ω) → L 4 (Ω), that is

uL 4  C u1 .

(3.1) h

h

We assume that finite element spaces have the following properties. The discrete spaces X , Q satisfy the usual approximation theoretic conditions and the inf–sup condition or Babuška–Brezzi condition i.e. there is a constant β independent of the mesh size h such that

inf

sup

qh ∈ Q h vh ∈ X h

(qh , ∇ · vh )  β > 0. ∇ vh qh 

(3.2)

For examples of such compatible spaces see e.g., [9,6]. Definition 3.1. Let V and V h denote respectively the divergence free subspaces of X and X h :





V := v ∈ X: (q, ∇ · v) = 0, ∀q ∈ Q ,









V h := vh ∈ X h : qh , ∇ · vh = 0, ∀qh ∈ Q h . Although typically V h  V , it is known that under the discrete inf–sup condition (3.2), functions in V are well approximated by ones in V h [6]. We consider X h and W k to be spaces of continuous piecewise polynomials of degree r and Q h is the space of continuous piecewise polynomials of degree r − 1. We also make the standard assumptions that the spaces X h , Q h and W k satisfy the following approximation properties for a given integer 1  s  r:

inf

vh ∈ X h , qh ∈ Q h

        u − vh + h ∇ u − vh + h p − qh  Ch s+1 us+1 +  p s ,





inf T − S k  ks+1  T s+1

(3.3) (3.4)

S k ∈W k

for (u, p , T ) ∈ ( X ∩ H s+1 (Ω), Q ∩ H s (Ω), W ∩ H s+1 (Ω)). We also use the fact that L 2 orthogonal projections of L H and M K satisfy

 G − P μ G   C μs | G |s ,

μ = H, K , 1  s  r

(3.5)

for G ∈ ( L (Ω) ∩ H (Ω)). We define the following weighted norms. 2

s

Definition 3.2. For u ∈ X , T ∈ W , the weighted norms of functions u : Ω f → R, T : Ω → R are defined by

2 ua2,b,α1 = au2 + b∇ u2 + α1 ( I − P H )∇ u , 2 = a T 2 + b∇ T 2 + α2 ( I − P K )∇ T  T 2 a,b,α2

where a, b > 0 are constants and α1 , α2 are stabilizing parameters. From now on, we denote min(κ f , κs ) as κmin and max(κ f , κs ) as

κmax for the sake of simplicity.

A. Çıbık, S. Kaya / J. Math. Anal. Appl. 381 (2011) 469–484

473

Lemma 3.1. The bilinear forms A 0 (·,·), A 1 (·,·) are continuous and coercive with respect to corresponding weighted norms. That is, for u, v ∈ X , T , S ∈ W , we have

A 0 (u, v)  u1,Pr,α1 v1,Pr,α1 , A 0 (u, u)  u21,Pr,α1 , A 1 ( T , S )   T 1,κmax ,α2  S 1,κmax ,α2 , A 1 ( T , T )   T 21,κmin ,α2 . Proof. Clearly, the bilinear forms (2.1) and (2.2) are continuous and coercive. The results follow from Definition 3.2.

2

Inequalities which are used frequently are Young’s inequality

ab 

t p

ap +

t −q/ p q

bq ,

a, b, p , q, t ∈ R,

1 p

+

1 q

= 1, p , q ∈ (1, ∞), t > 0,

(3.6)

and Poincaré’s inequality in X

v  C ∇ v

(3.7)

for all v ∈ X with C = C (Ω). Throughout this paper, the constant C is generic constant which depends on the domain Ω and independent from h, k, H , K , α1 and α2 unless stated otherwise. We remark on the convective terms defined by (2.4)–(2.5). In the continuous case, the standard form of the convective term and skew-symmetric form of trilinear form are identical if ∇ · u = 0 and if u vanishes on the boundary. Since standard convective terms are not divergence free on the finite element spaces, we use the modified ones [6]. Skew-symmetric trilinear forms (2.4)–(2.5) satisfy the following estimations. We denote by C 1 and C 2 finite positive constants with

c 0 (u, v, w)  C 1 ∇ u∇ v∇ w,

(3.8)

c 1 (u, T , S )  C 2 ∇ u∇ T ∇ S 

(3.9)

for all u, v, w ∈ X and T , S ∈ W . These estimates are well known and can be derived by applying Hölder’s inequality and Sobolev embedding theorems [6]. Remark 3.1. For u, v ∈ X and T ∈ W , we have c 0 (u, v, v) = 0 and c 1 (u, T , T ) = 0. Now, we also define the finite constant N h which used throughout the paper frequently:



 













N h = sup c 0 uh , vh , wh : ∇ vh = ∇ uh = ∇ wh = 1, uh , vh , wh ∈ V h . 4. Existence and uniqueness results of discrete problem Throughout this section, we consider the existence, uniqueness and stability properties of the discrete projection-based natural convection problem. These results without extra stabilization terms of continuous natural convection problem have been established in [3]. Using Lemma 3.1, similar results for continuous problem with stabilization can be established in the same way. For completeness, we only state and prove the existence and uniqueness of the discrete problem. Theorem 4.1 (Existence). The problem (2.12)–(2.13) has at least one solution. Proof. The proof consists of applying Lax–Milgram Theorem and Leray–Schauder Principle. Lax–Milgram Theorem guarantees the existence and uniqueness of T k in the solution of (2.13). Note that the approximate temperature T k depends on the velocity field uh . Thus we may define a mapping F hk : V h → W k by F hk (uh ) = T k . Now, we show that there is at least one uh ∈ V h satisfying











A 0 uh , vh + c 0 uh , uh , vh = Pr Ra d T k , vh



for all vh ∈ V h . From Lemma 3.1, A 0 (uh , vh ) is a continuous elliptic bilinear form on V h × V h and

        

−c 0 uh , uh , vh + Pr Ra d F hk uh  C ∇ uh 2 + Pr Ra F hk uh ∇ vh

(4.1)

474

A. Çıbık, S. Kaya / J. Math. Anal. Appl. 381 (2011) 469–484

for all vh ∈ V h . Thus, we may define a mapping G h : V h → V h by





















A 0 G h uh , vh = −c 0 uh , uh , vh + Pr Ra d F hk uh , vh . Note that uh is a solution of (4.1) if it is a solution of





G h uh = uh .

(4.2)

Thus, it suffices to show that there exists at least one solution to the fixed point problem (4.2). Leray–Schauder Principle guarantees the existence of a fixed point under two conditions: (i) G h should be completely continuous, (ii) there exists θ > 0 such that for every λ ∈ [0, 1] and vh ∈ V h with

  λG h vh = vh ,

(4.3)

vh should satisfy ∇ vh   θ . Since V h is finite-dimensional, G h is continuous and compact and thus completely continuous. This proves part (i). To prove the second condition, we consider only λ ∈ (0, 1] with λG h (vh ) = vh . Then, we have

        λ−1 A 0 vh , vh = −c 0 vh , vh , vh + Pr Ra d F hk vh , vh and

2 2   −1 λ−1 Pr ∇ vh + λ−1 α1 ( I − P H )∇ vh  Pr Ra ∇ F hk vh ∇ vh  Pr Ra κmin γ −1 ∇ vh . Hence

h ∇ v  λ Ra κ −1 γ −1 min

which completes the proof.

2

Before considering the uniqueness issue, we present some stability results. Lemma 4.1 (Stability of the velocity, temperature and pressure). The finite element approximation of (2.12)–(2.13) is stable in the following sense: −1 κmin ∇ T k 2 + 2α2 ( I − P K )∇ T k 2  κmin γ 2−1 , h 2 h 2 2 (ii) Pr ∇ u  + 2α1 ( I − P H )∇ u   Pr Ra  T k 2−1 , −2 (iii) Pr ∇ uh 2 + 2α1 ( I − P H )∇ uh 2  Pr Ra 2 κmin γ 2−1 , √ −1 h − 1 −1 (iv)  p   C β κmin γ −1 (Pr Ra + Pr α1 Ra + Ra 2 N h κmin γ −1 ).

(i)

Proof. To prove (i), we set S k = T k in (2.13) and apply the Young’s inequality. For (ii), we set uh = vh in (2.12) and use a similar argument as in (i). Combination of the parts (i) and (ii) gives (iii). To prove part (iv), consider Eq. (2.12) in X h :

















p h , ∇ · vh = A 0 uh , vh + c 0 uh , uh , vh − Pr Ra d T k , vh .

Cauchy–Schwarz inequality and (3.8) yield































2 p h , ∇ · vh  Pr ∇ uh ∇ vh + α1 ( I − P H )∇ uh ( I − P H )∇ vh + N h ∇ uh ∇ vh + Pr Ra T k −1 ∇ vh .

Making use of the stability bounds for the velocity and temperature gives

( ph , ∇ · vh ) −1 γ −1 +  Pr Ra κmin ∇ vh 



Pr α1 2

−1 −2 −1 γ −1 + N h Ra 2 κmin γ 2−1 + Pr Ra κmin γ −1 . Ra κmin

Taking supremum over vh ∈ X h and using the inf–sup condition (3.2) yield the desired result.

2

Corollary 4.1. Existence and uniqueness of p h is guaranteed by part (iv) of Lemma 4.1 and the inf–sup condition (3.2) [6]. We are now in a position to prove the global uniqueness condition of the discrete solution, which is the same as with the continuous case in [3]. First, by using the solution operator F hk in Theorem 4.1, we define the following constant:

 M hk = sup

d( F hk (uh ) − F hk (vh ), uh − vh )

∇(uh − vh )2

 ,u =

v , u ,v ∈ V . h

h

h

h

h

(4.4)

A. Çıbık, S. Kaya / J. Math. Anal. Appl. 381 (2011) 469–484

475

Theorem 4.2. Suppose N h ∇ uh  + Pr Ra M hk < Pr . Then, uh and F hk (uh ) = T k are unique solutions. Proof. Let uh , wh ∈ V h and uh = wh be two solutions. Writing Eq. (2.12) for uh and wh , and subtracting them give























A 0 uh − wh , vh = c 0 wh , wh , vh − c 0 uh , uh , vh + Pr Ra d F hk uh − F hk wh , vh



(4.5)

for all v ∈ V . Setting v = u − w in (4.5), using Cauchy–Schwarz inequality, adding and subtracting terms, using Remark 3.1 and (4.4) lead us to h

h

h



h

h



















2 2 2 2 Pr ∇ uh − wh + α1 ( I − P H )∇ uh − wh  N h ∇ uh − wh ∇ uh + Pr Ra ∇ uh − wh M hk .

So,









 









2 2 P r − N h ∇ uh + Pr Ra M hk ∇ uh − wh + α1 ( I − P H )∇ uh − wh  0.

Since ( N h ∇ uh  + Pr Ra M hk ) < Pr , we have a contradiction. Therefore, uh = vh .

2

Remark 4.1. If one uses the results of Lemma 4.1, global uniqueness condition, N h ∇ uh  + Pr Ra M hk < Pr can be reformu−1 −1 γ −1 ( N h + Pr C 2 κmin ) < Pr in terms problem data. lated as Ra κmin Furthermore, global uniqueness condition of the discrete problem ensures uh to be a fixed point of a contractive map in V h [17]. 5. A priori error estimation This section states a priori error estimation for the velocity and temperature. Before giving the main theorem, we define so-called modified Stokes projection operators. Lemma 3.1, hence Lax–Milgram Theorem guarantees the existence of such projection operators for both velocity and temperature. When we split the errors into approximation terms and a finite element remainder for u and T , the use of such operators simplifies the approximation terms and so the error estimations. We first state the stability of these projections and give the related error bounds. Definition 5.1 (Modified Stokes projections for the velocity and temperature). The operator of the modified Stokes projection for ˜ , p˜ ) where the velocity and pressure, P S , is defined by; P S : ( X , Q ) → ( X h , Q h ), P S (u, p ) = (u









˜ , vh + b vh , p − p˜ = 0, A0 u − u 



˜ , qh = 0 b u−u for all (vh , qh ) ∈ ( X h , Q h ). In the discretely divergence free space V h and in the pressure space Q h , this definition reduces to









˜ , vh + b vh , p − qh = 0 A0 u − u

(5.1)

for all v ∈ V . The modified Stokes projection operator for the temperature, P T , is defined by P T : W → W , P T ( T ) = T˜ where h

h



k



A 1 T − T˜ , S k = 0

(5.2)

for all S ∈ W . k

k

Lemma 5.1 (Stability of modified Stokes projections). The modified Stokes projections defined by (5.1) and (5.2) are stable in the following sense:



2





2



2 

˜ 2 + α1 ( I − P H )∇ u˜  C Pr ∇ u2 + α1 ( I − P H )∇ u + Pr −1 p − qh , Pr ∇ u

(5.3)

2

2 2 κ κmin ∇ T˜ 2 + α2 ( I − P K )∇ T˜  C max ∇ T 2 + α2 ( I − P K )∇ T .

(5.4)

κmin

˜ in (5.1) and use Cauchy–Schwarz inequality: Proof. For the proof of (5.3), first set vh = u



2











˜ 2 + α1 ( I − P H )∇ u˜  Pr ∇ u∇ u˜  + α1 ( I − P H )∇ u ( I − P H )∇ u˜ + p − qh ∇ · u˜ . Pr ∇ u Young’s inequality and combining terms give the result. The stability of the modified Stokes projection of temperature is established by writing S k = T˜ in (5.2) and using similar arguments as in the first part. 2 The next lemma states the error in those projection operators.

476

A. Çıbık, S. Kaya / J. Math. Anal. Appl. 381 (2011) 469–484

˜ , T˜ ) exists uniquely in Lemma 5.2 (Error in modified Stokes projections). Suppose the discrete inf–sup condition (3.2) holds. Then (u ( X h , Q h , W k ) and satisfies





2

2

˜ ) + α1 ( I − P H )∇(u − u˜ ) Pr ∇(u − u C















2 2 2 ˆ ) + α1 inf ( I − P H )∇(u − uˆ ) + Pr −1 inf p − qh , inf ∇(u − u

ˆ ∈Xh u

ˆ ∈Xh u



2

qh ∈ Q h

2

(5.5)

κmin ∇( T − T˜ ) + α2 ( I − P K )∇( T − T˜ )

 2 2  −1 2  C κmax κmin inf ∇( T − Tˆ ) + α2 inf ( I − P K )∇( T − Tˆ ) . Tˆ ∈ W k

˜ and decompose the error e = Proof. To prove (5.5), let e = u − u approximation of u in V h . Thus (5.1) reads as





(5.6)

Tˆ ∈ W k



Pr ∇φ h , ∇ vh + α1 ( I − P H )∇φ h , ( I − P H )∇ vh

ˆ φ = u˜ − u. ˆ Here uˆ is the η − φ h , where η = u − u,



      = Pr ∇ η, ∇ vh + α1 ( I − P H )∇ η, ( I − P H )∇ vh + p − qh , ∇ · vh .

(5.7)

Setting vh = φ in (5.7) and applying Cauchy–Schwarz and Young’s inequalities direct us to h



∇φ h 2 + α1 ( I − P H )∇φ h 2  C Pr ∇ η2 + α1 ( I − P H )∇ η 2 + Pr −1 p − qh 2 .

Pr 2

2

2

(5.8)

ˆ is an approximation of u in V h , we can take infimum over V h in (5.8). Recall that under the discrete inf–sup Since u condition (3.2) and ∇ · u = 0, the infimum can be replaced by X h [6]. The stated error estimate now follows from the triangle inequality. To prove (5.6), define T − T˜ = e˜ = ( T − Tˆ ) − ( T˜ − Tˆ ) = χ − ξ k where Tˆ approximates T in W k . As in the first part, if one sets S k = ξ k and uses Cauchy–Schwarz and Young’s inequalities, the following estimation is obtained



2



2



2

−1 2 κmin ∇ξ k + α2 ( I − P K )∇ξ k  C κmax κmin ∇ χ 2 + α2 ( I − P K )∇ χ .

Taking infimum over W k and applying the triangle inequality complete the proof.

2

We now give the our main theorems. Since the equations are coupled in (2.12)–(2.13), the error estimations are also coupled. Now we first state the error estimation for T − T k in terms of the error in u − uh . Theorem 5.1. The error for T − T k satisfies











2 2 κmin ∇ T − T k + α2 ( I − P K )∇ T − T k

 2   2  2 −1 −3  C κmin inf ∇ u2 ∇( T − T˜ ) + α2 ( I − P K )∇ T + C 22 κmin γ 2−1 ∇ u − uh . T˜ ∈ W k

Proof. Making use of (2.7) and (2.13) gives the error equation:















A 1 e˜ , S k + c 1 u, T , S k − c 1 uh , T k , S k = α2 ( I − P K )∇ T , ( I − P K )∇ S k



(5.9)

for all S k ∈ W k where e˜ = T − T k . Decompose the error as an approximation terms and a finite element remainder: e˜ = ( T − T˜ ) − ( T k − T˜ ) = χ − ξ k . Here, T˜ denotes the modified Stokes projection of T defined by (5.4). Now, set S k = ξ k into the error equation (5.9). With a rearrangement of terms, we obtain

 k k       

A 1 ξ , ξ  c 1 u, T , ξ k − c 1 uh , T k , ξ k + α2 ( I − P K )∇ T , ( I − P K )∇ξ k .

(5.10)

Note that A 1 (χ , ξ k ) = 0 due to the definition of the modified Stokes projection. Now, let us bound each term on the righthand side of (5.10):

       

c 1 u, T , ξ k − c 1 uh , T k , ξ k = c 1 u, χ , ξ k − c 1 u − uh , T k , ξ k

2  2 κmin k 2 −1 −1  ∇ξ ∇ u − uh ∇ T k , ∇ u2 ∇ χ 2 + + C 22 κmin  C κmin 4

 

α2 ( I − P K )∇ T , ( I − P K )∇ξ k  α2 ( I − P K )∇ T 2 + α2 ( I − P K )∇ξ k 2 . 2

2

A. Çıbık, S. Kaya / J. Math. Anal. Appl. 381 (2011) 469–484

477

Thus, bounding the terms as shown above for (5.10) results in

κmin 2

2 α2   ( I − P K )∇ξ k 2  C κ −1 ∇ u2 ∇ χ 2 + C 2 κ −1 ∇ u − uh 2 ∇ T k 2 + α2 ( I − P K )∇ T 2 . ∇ξ k + 2 min min 2

2

Combination of terms and application of the triangle inequality yield the stated error estimation.

2

The error estimation for the velocity is proved next. This error estimation uses Theorem 5.1. −1 −3 γ −1 ( N h + 32 C 22 Pr Ra κmin γ 3−1 ) < Pr , the error satisfies Theorem 5.2. Under the condition Ra κmin











2 2 Pr ∇ u − uh + α1 ( I − P H )∇ u − uh

  2  2 2  C M 1 inf ∇(u − u˜ ) + Pr −1 α1 inf ( I − P H )∇(u − u˜ ) + Pr −2 inf p − qh 

˜ ∈Xh u

−2 + M 2 κmin

˜ ∈Xh u

2 2  −1 inf ∇( T − T˜ ) + κmin α2 inf ( I − P K )∇( T − T˜ )

T˜ ∈ W k

qh ∈ Q h

T˜ ∈ W k

2 2  −1 α2 ( I − P K )∇ T , + α1 ( I − P H )∇ u + Pr Ra 2 κmin where C 2 is as in (3.9) and M 1 and M 2 are also constants which are defined below explicitly:





−2 −4 M 1 = C Pr −1 κmin Ra 2 γ 2−1 + Pr Ra 2 κmin γ 2−1 ,





κ2 M 2 = C Pr Ra max γ 2−1 . 6 κmin 4

Proof. The use of (2.7) and (2.12) results with the error equation:















A 0 e, vh + c 0 u, u, vh − c 0 uh , uh , vh + b vh , p − qh



    = Pr Ra d e˜ , vh + α1 ( I − P H )∇ u, ( I − P H )∇ vh

(5.11)

˜ ), φ h = (uh − u˜ ) for all (vh , qh ) ∈ ( V h , Q h ) where e = u − uh and e˜ = T − T k . Split the errors as e = η − φ h where η = (u − u k k k ˜ ˜ ˜ ˜ and T denote the modified Stokes projections of u and T , and e˜ = χ − ξ where χ = ( T − T ), ξ = ( T − T ). Note that u respectively. Now, writing vh = φ h in (5.11) yields





A0 φh , φh = A0















η, φ h + b φ h , p − qh + c0 u, u, φ h − c0 uh , uh , φ h



    + α1 ( I − P H )∇ u, ( I − P H )∇φ h + Pr Ra d e˜ , φ h .

(5.12)

Note that, A 0 (η, φ ) + b(φ , p − q ) = 0 by the definition of the modified Stokes projection. To bound the terms on the right-hand side of (5.12), we first consider the nonlinear terms. Adding, subtracting terms and observing the skew-symmetry of convective term yield h

h

h

         

c 0 u, u, φ h − c 0 uh , uh , φ h = c 0 u, η, φ h + c 0 η, uh , φ h − c 0 φ h , uh , φ h . Cauchy–Schwarz, (3.8) and Young’s inequalities give

 

c 0 u, η, φ h  C Pr −1 ∇ u2 ∇ η2 + Pr ∇φ h 2 , 6

 

c 0 η, uh , φ h  C Pr −1 ∇ η2 ∇ uh 2 + Pr ∇φ h 2 ,

 h h h 



c 0 φ , u , φ  N h ∇φ h 2 ∇ uh .

6

Similarly, consistency term and the last term on the right-hand side of (5.12) are bounded with

 

α1 ( I − P H )∇ u, ( I − P H )∇φ h  α1 ( I − P H )∇ u 2 + α1 ( I − P H )∇φ h 2 2

2

and

 

Pr Ra d T − T k , φ h  3 Pr Ra 2 T − T k 2 + Pr ∇φ h 2 . −1 2

6

h

Combining all the terms involving φ on the left-hand side gives

478

A. Çıbık, S. Kaya / J. Math. Anal. Appl. 381 (2011) 469–484



h h 2 α1 Pr ( I − P H )∇φ h 2 ∇φ + − Nh ∇ u 2

2

2  2  3 2   C Pr −1 ∇ η2 ∇ u2 + ∇ uh + α1 ( I − P H )∇ u + Pr Ra 2 T − T k −1 . 

(5.13)

2

Clearly, the next step we should follow is to find a bound for the term  T − T k 2−1 . In order to do that, we write u − uh =

η − φ h in the statement of Theorem 5.1 and plug in estimation in (5.13). Rearranging the terms yields

h 3 2 α1 Pr 2 2 −4 2 ( I − P H )∇φ h 2 − N h ∇ u − Pr Ra C 2 κmin γ −1 ∇φ h + 2

2

2

2  2 −4  C Pr ∇ η ∇ u + ∇ uh + α1 ( I − P H )∇ u + Pr Ra 2 κmin γ 2−1 ∇ η2 2  −2 −1 + Pr Ra 2 κmin ∇ u2 ∇ χ 2 + Pr Ra 2 κmin α2 ( I − P K )∇ T . 

−1

2



2

(5.14)

Let us consider the coefficient of the term ∇φ  . Making use of the uniqueness bound and the assumption of the theorem, we have h 2

Pr 2

<

Pr 2

−1 − N h Ra κmin γ −1 −

3 2

−4 γ 2−1 . Pr Ra 2 C 22 κmin

(5.15)

Plugging (5.15) into (5.14) and writing the stability bounds for the terms we have



2



2



−2 −4 Pr ∇φ h + α1 ( I − P H )∇φ h  C Pr −1 κmin Ra 2 γ 2−1 ∇ η 2 + Pr Ra 2 κmin γ 2−1 ∇ η2

2 −4 + Pr Ra 4 κmin γ 2−1 ∇ χ 2 + α1 ( I − P H )∇ u 2  −1 + Pr Ra 2 κmin α2 ( I − P K )∇ T .

Substituting the error bounds of Lemma 5.2 into (5.16) and applying the triangle inequality complete the proof.

(5.16)

2

One might see also that the addition of the extra term in (2.12)–(2.13) does not degrade the order of convergence. To see this, we give the following remark. Remark 5.1. If we assume the regularity assumptions, (u, p , T ) ∈ ( X ∩ H s+1 (Ω), Q ∩ H s (Ω), W ∩ H s+1 (Ω)) and the use of the estimations (3.3), (3.4) and (3.5) yield











2 2 Pr ∇ u − u h + α1 ( I − P H )∇ u − u h

     M 1 h2s |u |2s+1 1 + P r −1 α1 + P r −2 h2s | p |2s  −1 2s 2  −1   k | T |s+1 κmin + α2 + α1 H 2s |u |2s+1 + α2 K 2s | T |2s+1 . + M 2 κmin

(5.17)

Here h, k are given and by equilibrating the orders of convergence, appropriate values for the mesh scales H , K and parameters α1 , α2 are chosen. That is, the error is optimal for α1 H 2s = h2s and α2 K 2s = k2s . For instance, let us consider the case for s = 2 and use Taylor–Hood finite element pairs, satisfying the inf–sup condition (3.2), which are given below explicitly along with the choices of L H = ∇ X H and M K = ∇ W K :





¯ : v| ∈ P 2 (), ∀ ∈ F h , X h = v ∈ C 0 (Ω) 







¯ : S | ∈ P 2 (), ∀ ∈ G k , W k = S ∈ C 0 (Ω) ¯ : v| ∈ P 1 (), ∀ ∈ F h , Q h = v ∈ C 0 (Ω) 





L H = l H ∈ L 2 (Ω): l H  ∈ P 1 (), ∀ ∈ F H ,







M K = m K ∈ L 2 (Ω): m K  ∈ P 1 (), ∀ ∈ G K . If we make the same assumptions as in Theorem 5.2 and consider (5.17), one can imply that along with the choices of

(α1 , H ) = (h2 , h1/2 ) the error

    ∇ u − uh = O h2 + k2

is optimal for the velocity. If we plug the results obtained for velocity into Theorem 5.1 and carry out similar computations, we have

    ∇ T − T k = O h 2 + k 2 . Similarly, for the choices of (α2 , K ) = (k2 , k1/2 ) we have the optimal error for the temperature.

A. Çıbık, S. Kaya / J. Math. Anal. Appl. 381 (2011) 469–484

479

6. Error estimation for pressure This section deals with the estimation of the error for the discrete pressure in the L 2 norm. Theorem 6.1 (Error estimate for pressure). Suppose that the assumptions of Theorem 5.2 hold. Then the error p − p h satisfies

        p − ph ≤ qC Pr + ∇ u ∇ u − uh + ∇ u − uh 2 + α1 ( I − P H )∇ u − uh inf p − qh qh ∈ Q h

   + Pr Ra T − T k −1 + α1 ( I − P H )∇ u .

Proof. To prove this, we consider (5.11). Let p − p h = ( p − p˜ ) − ( p h − p˜ ), where p˜ is an approximation of the pressure in Q h . (5.11) reads as









 





b vh , p h − p˜ = A 0 e, vh + c 0 u, u, vh − c 0 uh , uh , vh

  − α1 ( I − P H )∇ u, ( I − P H )∇ vh .



    + b vh , p − p˜ − Pr Ra d e˜ , vh

We first consider here the nonlinear terms. Adding and subtracting terms and using (3.8) yield

   

     

c 0 u, u, vh − c 0 uh , uh , vh = −c 0 e, e, vh + c 0 e, u, vh + c 0 u, e, vh

   C ∇ e + ∇ u ∇ e ∇ vh . Bounds for the other terms are obtained in a similar manner as in the estimation of the error ∇(u − uh ). Hence

 h h   



b v , p − p˜  C ∇ vh Pr ∇ e + α1 ( I − P H )∇ e + ∇ e + ∇ u ∇ e  +  p − p˜  + Pr Ra T − T k −1 + α1 ( I − P H )∇ u .

(6.1)

Notice that (3.2) implies











p h − p˜ , ∇ · vh  β ∇ vh p h − p˜

and using this relation yields h h p − ph   p − p˜  + p˜ − ph   p − p˜  + β −1 |b(v , p − p˜ )| . ∇ vh 

Substituting (6.1) into (6.2) taking infimum over Q h give us the desired result.

(6.2)

2

Remark 6.1. Making use of Taylor–Hood elements as in Remark 5.1 with the choices (α1 , H ) = (h2 , h1/2 ) and (α2 , K ) = (k2 , k1/2 ) and using the approximation results (3.3)–(3.4) for the velocity and temperature errors, we have

  p − p h = O h 2 + k2 which is the optimal error. 7. Numerical studies In this section, numerical studies are given in order to show the effectiveness of the method and validate the obtained theoretical results. The projection-based stabilization method for steady natural convection problem has been assessed on two numerical examples in two dimensions. The first example is a well-known test case for the natural convection codes which is called buoyancy-driven cavity problem. For the other test case, we consider a known particular analytical solution in order to check the error rates. All computations are carried out by using the software FreeFem++ [10]. In both examples, we use conforming Taylor–Hood finite element pairs. It is well known that these pairs fullfill the inf–sup condition (3.2) (see [9]). Finite element spaces are given in Remarks 5.1 and 6.1 with the algorithmic choices for the size of the meshes and the parameters: H ∼ h 1/2 and K ∼ k1/2 , α1 = h2 , α2 = k2 . Since we solve the problem on the same mesh, we let h = k and H = K . To handle the nonlinearity of the system, the Newton method of [9] is used. The algorithm consists of starting with an initial guess (u(0) , T (0) ) and then generate the sequence of iterates (u(m) ∈ X h , p (m) ∈ Q h and T (m) ∈ W k ) for m  1 by solving the sequence of linear systems

480

A. Çıbık, S. Kaya / J. Math. Anal. Appl. 381 (2011) 469–484

Fig. 1. The physical domain with its boundary conditions. Table 1 Comparison of maximum vertical velocity at y = 0.5 with mesh size used in computation. Ra

GFEM

Present study

Ref. [5]

Ref. [21]

Ref. [20]

Ref. [26]

104

16.41(11 × 11)

19.91(11 × 11)

19.51(41 × 41)

19.63(71 × 71)

19.90(71 × 71)

105

51.22(21 × 21)

70.60(21 × 21)

68.22(81 × 81)

68.85(71 × 71)

70.00(71 × 71)

70.63(101 × 101)

106

201.20(32 × 32)

228.12(32 × 32)

216.75(81 × 81)

221.60(71 × 71)

228.00(71 × 71)

227.11(101 × 101)















Pr a0 u(m) , vh + c 0 u(m−1) , u(m) , vh + c 0 u(m) , u(m−1) , vh + b vh , p (m)

19.79(101 × 101)



        = Pr Ra d T (m) , vh + c 0 u(m−1) , u(m−1) , vh − α1 ( I − P H )∇ u(m−1) , ( I − P H )∇ vh b u(m) , qh = 0,       a1 T (m) , S k + c 1 u(m) , T (m−1) , S k + c 1 u(m−1) , T (m) , S k       = γ , S k + c 1 u(m−1) , T (m−1) , S k − α2 ( I − P H )∇ T (m−1) , ( I − P H )∇ S k for all (vh , qh , S k ) ∈ ( X h , Q h , W k ). This scheme is known to be locally convergent if at least either or both T and u · n are specified at every point of the boundary. 7.1. Buoyancy-driven cavity problem The problem of buoyancy-driven cavity is used as a suitable benchmark for testing the natural convection codes in the literature. The simplicity of geometry and clear boundary conditions make this problem attractive. The domain consists of a square cavity with differentially heated vertical walls where right and left walls are kept at T C and T H , respectively, with T H > T C . The remaining walls are insulated and there is no heat transfer through them. The boundary conditions are no-slip boundary conditions for the velocity at all four walls (u = 0) and Dirichlet boundary conditions for the temperature at vertical walls. As the horizontal walls are adiabatic, we employ ∂∂ Tn = 0 here. Fig. 1 shows the physical domain of the buoyancy-driven cavity flow problem. In the test case, we take κ = 1, γ = 0, T C = 0 and T H = 1. While we consider the air as the cavity filling fluid in our model, we take the fixed value Pr = 0.71. We have performed our computations for Rayleigh number varying from 103 to 106 . The performance of the projection-based stabilization is compared with the famous benchmark solutions of de Vahl Davis [5] and some other authors such as Massarotti et al. [21], Manzari [20], and the more recent study of Wan et al. [26]. From these benchmark solutions [5] used second-order central approximations to solve natural convection problem in a square cavity. Ref. [21] developed a semi-implicit form of the characteristic-based split scheme and [20] employed an explicit finite element algorithm. Recently, [26] used discrete singular convolution for the solution of the problem. We also include the results for the classical Galerkin Finite Element Method (GFEM) where we keep the same mesh sizes for the proposed method and GFEM. Numerical simulations are obtained for three uniform grids of 11 × 11, 21 × 21 and 32 × 32. We start our illustrations by giving peak values of vertical velocity at y = 0.5 and horizontal velocity at x = 0.5. Tables 1 and 2 summarize the maximum vertical velocity values at mid-height and at mid-width for different Rayleigh numbers. For quantitative assessment, we also include those velocity values obtained by [5,21,20,26]. As can be observed, the results of our computations are in an excellent agreement with the benchmark data even at coarser grid. We also see that as the Rayleigh number increases, GFEM yields results which are not so close to the benchmark solutions. We also present the vertical velocity distribution at the mid-height and horizontal velocity distribution at the mid-width in Fig. 2, respectively,

A. Çıbık, S. Kaya / J. Math. Anal. Appl. 381 (2011) 469–484

481

Table 2 Comparison of maximum horizontal velocity at x = 0.5 with mesh size used in computation. Ra

GFEM

Present study

Ref. [5]

Ref. [20]

Ref. [26]

104

15.70(11 × 11)

15.90(11 × 11)

16.18(41 × 41)

16.10(71 × 71)

16.10(101 × 101)

105

41.00(21 × 21)

33.51(21 × 21)

34.81(81 × 81)

34.00(71 × 71)

34.00(101 × 101)

106

80.25(32 × 32)

65.52(32 × 32)

65.33(81 × 81)

65.40(71 × 71)

65.40(101 × 101)

Fig. 2. Variation of vertical velocity at mid-height (left) and horizontal velocity at mid-width for varying Rayleigh numbers (right). Table 3 Comparison of average Nusselt number on the vertical boundary of the cavity at x = 0 with mesh size used in computation. Ra

GFEM

Present study

Ref. [5]

Ref. [20]

Ref. [21]

Ref. [26]

104

2.40(11 × 11)

2.15(11 × 11)

2.24(41 × 41)

2.08(71 × 71)

2.24(71 × 71)

2.25(101 × 101)

105

5.11(21 × 21)

4.35(21 × 21)

4.52(81 × 81)

4.30(71 × 71)

4.52(71 × 71)

4.59(101 × 101)

106

6.00(32 × 32)

8.83(32 × 32)

8.92(81 × 81)

8.74(71 × 71)

8.82(71 × 71)

8.97(101 × 101)

which are very popular graphical illustrations in the study of buoyancy-driven cavity type tests. These profiles are also comparable with the similar ones in [26]. It is obvious that as Rayleigh numbers increases, the differences in the profiles presented in Fig. 2 are getting larger. A very important property of the natural convection flows, especially for engineers, is the rate of heat transfer along the vertical walls of the cavity. The dimensionless parameter called Nusselt number stands for this quantity. The local Nusselt number can be calculated as

Nulocal = ±

∂T . ∂x

The negative sign means heat transfer at the hot wall and the positive sign means heat transfer at the cold wall. The local Nusselt number at the cavity hot wall is used for comparison with benchmark problems in the literature frequently. As in the velocity components case, we calculate the average Nusselt numbers with GFEM and our method. The benchmark data results are also included to compare the average Nusselt numbers values with the presented study. The results are given in Table 3. As we can understand from Table 3, there is a very good agreement with the benchmark solutions and the present study, which can still capture reasonable results for rather coarser grid. The plots of Fig. 3 show the variation of the Nusselt number along the hot wall and cold wall for different Rayleigh numbers. These profiles are also look reasonable when compared with those reported in [5,20,21,26]. Characters of the flow patterns for increasing Rayleigh numbers are seen very often in the study of natural convection problems. Diagrams showing the streamlines and temperature isolines are very popular among the convective heat transport illustrations. We present these patterns in Fig. 4. It is clear from the

482

A. Çıbık, S. Kaya / J. Math. Anal. Appl. 381 (2011) 469–484

Fig. 3. Variation of local Nusselt number at cavity hot wall (left) and cavity cold wall (right).

Fig. 4. Streamlines (upper left to right) and isotherms (lower left to right) for with Ra = 103 , 104 , 105 , 106 , respectively.

streamline patterns that, as Rayleigh number increases circular vortex at the cavity center begin to deform into an ellipse and then break up into two vortices tending to approach to the corners differentially heated sides of the cavity. So we can conclude that, the flow is swifter as the thermal convection is concentrated. Through the increase in Rayleigh number, parallel behavior of the temperature isolines is distorted and these lines seem to have a flat behavior in the central part of the region. Near the sides of the cavity, isolines tend to be vertical only. With Ra = 106 , the temperature slopes at the corners of the differentially heated sides are more immersed then the cases of lower Rayleigh number. We also note that these graphics are also perfectly comparable with the ones given in the investigations of [5,21,20,26].

A. Çıbık, S. Kaya / J. Math. Anal. Appl. 381 (2011) 469–484

483

Table 4 Total degree of freedoms, numerical errors and convergence rates for each variable.

u − uh 

Rate

∇(u − uh )

Rate

 p − ph 

Rate

∇( T − T k )

h = 1/4

374

0.0170

–

0.3712

–

0.3521

–

0.2922

Mesh

# of d.o.f.

Rate

h = 1/8

1318

0.0021

2.85

0.0905

2.02

0.0951

1.92

0.0767

1.95

h = 1/16

4934

2.434e-04

2.92

0.0222

2.01

0.0215

1.94

0.0187

2.02

h = 1/32

19078

2.722e-05

2.99

0.0054

2.01

0.0054

1.98

0.0042

2.10

7.2. Numerical convergence study An assessment of the convergence of the numerical simulation is presented in this subsection. We consider the problem (1.1) in the domain Ω = [−1, 1] × [−1, 1]. The forcing function γ and boundary values of the temperature are given so that the prescribed solution of the system is given by

u=



x2 − 1

1

p=

8

+

2 





 

y y 2 − 1 , − x2 − 1 x y 2 − 1



y 4 − y 6 + y 2 − 1 x8 +

1

6

4

2



1 2

2

2 

, 

y 6 − y 4 − y 2 + 1 x6 +



3



6 5

yx5 +

3 4







y 4 − y 6 + y 2 − 1 x4 + 4 y 3 − 8 y x3

y − y − y + 1 x + 10 y − 4 y x, 2

   1  3 1  5 3 5 3 3 1 1  2 T= 2 y − 3 y 5 + y x8 + 3 y − 2 y 3 − y x6 + y 3 − y 5 + y x4 + 3 y − 2 x3 x + 400 100 250 100 2 2 25

+

1  100



3 y 5 − 2 y 3 − y x2 +

1  50



3 y 4 − 12 y 2 + 8 x.

In (1.1), non-zero Neumann boundary condition for T on Γ B and Dirichlet boundary condition for u are chosen so that (u, p , T ) is the solution of the system. Note that, using the non-zero Neumann boundary condition for the variable T affect the stability bounds given in Lemma 4.1 and hence the main theorems. Although this replacement changes some terms and constants in the error analysis, it does not degrade the order of errors given in Remarks 5.1 and 6.1. We use the same settings as in Remarks 5.1 and 6.1 with Pr = 1, Ra = 100 and κ = 1. We compute the errors between exact solution and computed numerical solution for the variables u, p and T . Then, we compare error rates with the theoretical expectations. Table 4 presents the corresponding total degree of freedoms for u, T and p, errors and convergence rates for different mesh sizes. We first compute the errors for the coarsest mesh of h = 1/4 and then refine the mesh to obtain finer ones. The theory predicts the error rates in Table 4, O (h3 ) for the L 2 norm for u, and O (h2 ) for the L 2 norm for p and O (h2 ) in energy norm for the temperature. Note that the behavior of the error is exactly as anticipated by the theory. Thus we can conclude that the projection-based stabilization does not degrade the order of the errors for all variables. 8. Conclusion This paper studied the projection-based stabilization method for the steady-state natural convection equations. By means of this method, global stabilizations are added for both velocity and temperature variables and these effects are subtracted from the large scales. We established the rigorous finite element error analysis of the scheme for the velocity, temperature and pressure and proved that with the appropriate choices of mesh scales and the parameters, the optimal errors can be obtained. We examined performance and accuracy of the method and compared the results with other published data. The numerical results revealed excellent agreement with other published data and validation of theoretical results. References [1] R.A. Adams, Sobolev Spaces, Academic Press, New York, 1975. [2] J. Boland, W. Layton, An analysis of the finite element method for natural convection problems, Numer. Methods Partial Differential Equations 2 (1990) 115–126. [3] J. Boland, W. Layton, Error analysis for finite element methods for steady natural convection problems, Numer. Funct. Anal. Optim. 11 (1990) 449–483. [4] R. Codina, Comparison of some finite element methods for solving the diffusion–convection–reaction equation, Comput. Methods Appl. Mech. Engrg. 156 (1998) 185–210. [5] D. de Vahl Davis, Natural convection of air in a square cavity: A benchmark solution, Internat. J. Numer. Methods Fluids 3 (1983) 249–264. [6] V. Girault, P.A. Raviart, Finite Element Approximation of the Navier–Stokes Equations, Lecture Notes in Math., vol. 749, Springer-Verlag, Berlin, 1979. [7] P.M. Gresho, M. Lee, S.T. Chan, R.L. Sani, Solution of time dependent, incompressible Navier–Stokes and Boussinesq equations using the Galerkin finite element method, in: Lecture Notes in Math., vol. 771, Springer-Verlag, Berlin/Heidelberg/New York, 1980, pp. 203–222. [8] J.-L. Guermond, Stabilization of Galerkin approximations of transport equations by subgrid modeling, M2AN 33 (1999) 1293–1316. [9] M. Gunzburger, Finite Element Methods for Viscous Incompressible Flows: A Guide to Theory, Practice and Algorithms, Academic Press, Boston, 1989. [10] F. Hecht, FreeFem++ manual, http://www.freefem.org/ff++, 2009. [11] N. Heitmann, Subgrid stabilization of time-dependent convection dominated diffusive transport, J. Math. Anal. Appl. 331 (2007) 38–50.

484

A. Çıbık, S. Kaya / J. Math. Anal. Appl. 381 (2011) 469–484

[12] T.J.R. Hughes, Multiscale phenomena: Green’s functions, the Dirichlet-to-Neumann formulation, subgrid-scale models bubbles and the origin of stabilized methods, Comput. Methods Appl. Mech. Engrg. 127 (1995) 387–401. [13] V. John, S. Kaya, A finite element variational multiscale method for the Navier Stokes equations, SIAM J. Sci. Comput. 26 (2005) 1485–1503. [14] V. John, S. Kaya, W. Layton, A two-level variational multiscale method for convection–diffusion equations, Comput. Methods Appl. Mech. Engrg. 195 (2005) 4594–4603. [15] S. Kaya, B. Riviere, A two-grid stabilization method for solving the steady-state Navier–Stokes equations, Numer. Methods Partial Differential Equations 3 (2006) 728–743. [16] S.A. Korpela, D. Gözüm, C.B. Baxi, On the stability of the conduction regime or natural convection in a vertical slot, Internat. J. Heat Mass Transfer 16 (1973) 1683–1690. [17] W. Layton, Introduction to Finite Element Methods for Incompressible, Viscous Flows, SIAM Publ., 2008. [18] W.J. Layton, A connection between subgrid scale eddy viscosity and mixed methods, Appl. Math. Comput. 133 (2002) 147–157. [19] J. Löwe, G. Lube, A projection-based variational multiscale method for large-eddy simulation with application to non-isothermal free-convection problems, Technical report, NAM-Preprint, 2010. [20] M.T. Manzari, An explicit finite element algorithm for convective heat transfer problems, Internat. J. Numer. Methods Heat Fluid Flow 9 (1999) 860–877. [21] N. Massarotti, P. Nithiarasu, O.C. Zienkiewicz, Characteristic-Based-Split (CBS) algorithm for incompressible flow problems with heat transfer, Internat. J. Numer. Methods Heat Fluid Flow 8 (1998) 969–990. [22] M. Braack, E. Burman, V. John, G. Lube, Stabilized finite element methods for the generalized Oseen problem, Comput. Methods Appl. Mech. Engrg. 196 (2007) 853–866. [23] H. Melhem, Finite element approximation to heat transfer through combined solid and fluid media, PhD thesis, University of Pittsburgh, 1987. [24] P.H. Rabinowitz, Existence and nonuniqueness of rectangular solutions of the Benard problem, Arch. Ration. Mech. Anal. 29 (1968) 32–57. [25] H.-G. Roos, M. Stynes, L. Tobiska, Robust Numerical Methods for Singularly Perturbed Differential Equations, Springer Ser. Comput. Math., SpringerVerlag, Berlin, 2008. [26] D.C. Wan, B.S.V. Patnaik, G.W. Wei, A new benchmark quality solution for the buoyancy-driven cavity by discrete singular convolution, Numer. Heat Transfer, Part B 40 (2001) 199–228.