Robust globally divergence-free weak Galerkin finite element methods for natural convection problems

# Robust globally divergence-free weak Galerkin finite element methods for natural convection problems

[
###### Abstract

This paper proposes and analyzes a class of weak Galerkin (WG) finite element methods for stationary natural convection problems in two and three dimensions. We use piecewise polynomials of degrees and for the velocity, pressure, and temperature approximations in the interior of elements, respectively, and piecewise polynomials of degrees for the numerical traces of velocity, pressure and temperature on the interfaces of elements. The methods yield globally divergence-free velocity solutions. Well-posedness of the discrete scheme is established, optimal a priori error estimates are derived, and an unconditionally convergent iteration algorithm is presented. Numerical experiments confirm the theoretical results and show the robustness of the methods with respect to Rayleigh number.

Han Y H et. al]Yihui Han, Xiaoping Xie111Corresponding author.

• Robust globally divergence-free weak Galerkin finite element methods for natural convection problems

• [

• School of Mathematics, Sichuan University, Chengdu 610064, China

• AMS subject classifications: 52B10, 65D18, 68U05, 68U07

• Key words: natural convection, Weak Galerkin method, Globally divergence-free, error estimate, Rayleigh number.

## 1 Introduction

Let be a polygonal or polyhedral domain with a polygonal or polyhedral subdomain and , we consider the following stationary natural convection (or conduction-convection) problem: seek the velocity , the pressure , and the temperature such that

 ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩−PrΔu+∇⋅(u⊗u)+∇p−PrRajT=finΩf,∇⋅u=0inΩf,−κΔT+∇⋅(uT)=ginΩ,u≡0inΩs⋃∂Ωf,T=0on∂Ω, (1.1)

where is defined by for , is the vector of gravitational acceleration with when and when , , are the forcing functions, and , denote the Prandtl and Rayleigh numbers, respectively,.

The model problem (1.1), arising both in nature and in engineering applications, is a coupled system of fluid flow, governed by the incompressible Navier-Stokes equations, and heat transfer, governed by the energy equation. Due to its practical significance, the development of efficient numerical methods for natural convection has attracted a great many of research efforts; see, e.g. [1],[2],[3],[29],[23],[12],[15],[18],[20],[19],[21],[22],[25],[26],[28],[38],[39]. In [2, 3], error estimates for some finite element methods were derived in approximating stationary and non-stationary natural convection problems. [20, 19] applied Petrov-Galerkin least squares mixed finite element methods to discretize the problems. [25, 26] developed a nonconforming mixed element method and a Petrov-Galerkin least squares nonconforming mixed element method for the stationary problems. In [37], three kinds of decoupled two level finite element methods were presented. [38, 39] applied the variational multiscale method to solve the stationary and non-stationary problems.

In this paper, we consider a weak Galerkin (WG) finite element discretization of the model problem (1.1). The WG method was first proposed and analyzed to solve second-order elliptic problems [30, 31]. It is designed by using a weakly defined gradient operator over functions with discontinuity, and then allows the use of totally discontinuous functions in the finite element procedure. Similar to the hybridized discontinuous Galerkin (HDG) method [11], the WG method is of the property of local elimination of unknowns defined in the interior of elements. We note that in some special cases the WG method and the HDG method are equivalent (cf. [6, 7, 8]). In [6], a class of robust globally divergence-free weak Galerkin methods for Stokes equations were developed, and then were extended in [40] to solve incompressible quasi-Newtonian Stokes equations.We also refer to [9, 16, 13, 17, 24, 33, 32, 41, 10, 35, 34, 36] for some other developments and applications of the WG method.

This paper aims to propose a class of WG methods for the natural convection problems. The methods include as unknowns the velocity, pressure, and temperature variables both in the interior of elements and on the interfaces of elements. In the interior of elements, we use piecewise polynomials of degrees and for the velocity, pressure, and temperature approximations, respectively. On the interfaces of elements, we use piecewise polynomials of degrees for the numerical traces of velocity, pressure and temperature. The methods are shown to yield globally divergence-free velocity approximations.

The rest of the paper is organized as follows. Section 2 introduces the WG finite element scheme. Section 3 shows the existence and uniqueness of the discrete solution. Section 4 derives a priori error estimates. Section 5 discusses the local elimination property and the convergence of an iteration method for the WG scheme. Finally, Section 6 provides numerical examples to verify the theoretical results.

Throughout this paper, we use to denote , where the constant C is positive independent of mesh size and the , and Rayleigh number.

## 2 WG finite element scheme

### 2.1 Notation

For any bounded domain , let and denote the usual -order Sobolev spaces on D, and denote the norm and semi-norm on these spaces. We use to denote the inner product of , with . When , we set and . In particular, when , we use to replace . For integer , denotes the set of all polynomials on D with degree no more than . We also need the following spaces:

,

Let and be shape-regular simplicial decompositions of the subdomains and , respectively. Then is a shape-regular simplicial decomposition of . Let and be the sets of all edges (faces) of all elements in and , respectively, and set . For any , , we denote by and the diameters of and , respectively, and set . Let and be the outward unit normal vectors along the boundary and . We denote by and the piecewise-defined gradient and divergence with respect to . We also introduce the mesh-dependent inner products and mesh-dependent norms:

 ⟨u,v⟩∂Th:= ∑K∈Th⟨u,v⟩∂K,∥u∥0,∂Th:=⎛⎝∑K∈Th∥u∥20,∂K⎞⎠1/2 (u,v)Th:= ∑K∈Th(u,v)K,∥u∥0,Th:=⎛⎝∑K∈Th∥u∥20,K⎞⎠1/2.

### 2.2 Weak problem

We first introduce the space

 W:={v∈[H10(Ωf)]d:∇⋅v=0}

and the following bilinear and trilinear forms: for any , , and ,

 a(u,v) :=Pr(∇u,∇v),b(v,q):=(∇q,v), d(T,v) :=PrRa(jT,v),¯¯¯a(T,s):=κ(∇T,∇s), c(w;u,v) :=((w⋅∇)u,v),¯¯c(w;T,s):=((w⋅∇)T,s).

It is easy to see that, for ,

 c(u;u,v) =12(∇⋅(u⊗u),v)−12(∇⋅(u⊗v),u), ¯¯c(u;T,s) =12(∇⋅(uT),s)−12(∇⋅(us),T).

Then the variational problem of (1.1) reads as follows: seek such that

 {A(u;u,v)+b(v,p)−b(u,q)−d(T,v)=(f,v),∀(v,q)∈W×L20(Ωf),¯¯¯¯A(u;T,s)=(g,s),∀s∈H10(Ω), (2.1)

where

 A(u;u,v) :=a(u,v)+c(u;u,v), ¯¯¯¯A(u;T,s) =¯¯¯a(T,s)+¯¯c(u;T,s).
###### Theorem 2.1.

[3]For and , the weak problem (2.1) has at least one solution. In addition, it admits a unique solution if

 (Pr−1NRaκ−1+RaMκ−2)∥g∥−1+NPr−2∥f∥−1<1,

where

 N:=sup0≠w,u,v∈Wc(w;u,v)|w|1|u|1|v|1,M:=sup0≠w∈W,0≠T,s∈H10(Ω)¯¯c(w;T,s)|w|1|T|1|s|1.

In what follows, we assume that the solution is unique and, more precisely, there exists a fixed constant such that

 (Pr−1NRaκ−1+RaMκ−2)∥g∥−1+NPr−2∥f∥−1<1−δ.

### 2.3 Discrete weak operators

In order to design a WG finite element scheme for the problem (1.1), we introduce the discrete weak gradient operator and the discrete weak divergence operator as follows.

###### Definition 2.1.

For any and , the discrete weak gradient on is determined by the equation

 (∇w,r,Kv,τ)K=−(v0,∇⋅τ)K+⟨vb,τ⋅nK⟩∂K∀τ∈[Pr(K)]d.

Then we define the global discrete weak gradient operator by

.

For a vector , we define its discrete weak gradient by

 ∇w,rv:=(∇w,rv1,⋯,∇w,rvd)T.
###### Definition 2.2.

For any and , the discrete weak divergence is determined by the equation

 (∇w,r,K⋅v,τ)K=−(v0,∇τ)K+⟨vb⋅nK,τ⟩∂K∀τ∈Pr(K).

Then we define the global discrete weak divergence operator by

 ∇w,r⋅|K=∇w,r,K⋅,∀K∈Th.

For a tensor with for , we define its discrete weak divergence by

 ∇w,r⋅w=(∇w,r⋅w1,⋯,∇w,r⋅wd)T.

### 2.4 WG finite element scheme

For any and any integer , let and be the usual projection operators. We shall use to denote for vector spaces.

For any integer and , we introduce the following finite dimensional spaces:

 Vh ={vh={vh0,vhb}:vh0|K∈[Pk(K)]d,vhb|e∈[Pl(e)]d,∀K∈Th,∀e∈εh}, V0h ={vh={vh0,vhb}∈Vh:vhb|∂Ωf=0}, Qh ={qh={qh0,qhb}:qh0|K∈Pk−1(K),qhb|e∈Pk(e),∀K∈Th,∀e∈εh}, Q0h ={qh={qh0,qhb}∈Qh:qh0∈L20(Ωf)}, Sh ={sh={sh0,shb}:sh0|K∈Pk(K),shb|e∈Pl(e),∀K∈Th,∀e∈εh}, S0h ={sh={sh0,shb}∈Sh:shb|∂Ω=0}.

For any , , and , define the following bilinear and trilinear forms:

 ah(uh,vh) :=Pr(∇w,muh,∇w,mvh)+Pr⟨τ(Qbluh0−uhb),Qblvh0−vhb⟩∂Tfh, bh(vh,qh) :=(∇w,kqh,vh0), dh(Th,vh) :=PrRa(jTh0,vh0), ¯¯¯ah(Th,sh) :=κ(∇w,mTh,∇w,msh)+κ⟨τ(QblTh0−Thb),Qblsh0−shb⟩∂Th, ch(wh;uh,vh) :=12(∇w,k⋅(uh⊗wh),vh0)−12(∇w,k⋅(vh⊗wh),uh0), ¯¯ch(uh;Th,sh) :=12(∇w,k⋅(uhTh),sh0)−12(∇w,k⋅(uhsh),Th0).

It is easy to see that

 ch(wh;vh,vh)=0,¯¯ch(wh;sh,sh)=0. (2.2)

The WG finite element scheme for (1.1) is then given as follows: seek , , and such that

 {Ah(uh;uh,vh)+bh(vh,ph)−bh(uh,qh)−dh(Th,vh)=(f,vh0),∀(vh,qh)∈V0h×Q0h,¯¯¯¯Ah(uh;Th,sh)=(g,sh0),∀sh∈S0h, (2.3)

where

 Ah(wh;uh,vh) :=ah(uh,vh)+ch(wh;uh,vh), (2.4) ¯¯¯¯Ah(uh;Th,sh) :=¯¯¯ah(Th,sh)+¯¯ch(uh;Th,sh), (2.5)

, and m is an integer with .

###### Remark 2.1.

It’s easy to show that the scheme (2.3) yields globally divergence-free velocity approximation . In fact, let be any two adjacent elements with a common face , introduce a function with

 rhb|e={−(uh0⋅ne)|K1⋂e−(uh0⋅ne)|K2⋂e,∀e∈εfh/∂Ωf,0,∀e∈∂Ωf,

and set . Then, taking in (2.3) yields

 ∥∇h⋅uh0∥20+∑e∈εfh/∂Ωf∥(uh0⋅ne)|K1+(uh0⋅ne)|K2∥20,e=0.

This indicates and , i.e. the velocity approximation is globally divergence-free in a pointwise sense.

## 3 Well-posedness of the discrete scheme

### 3.1 Some basic results

For the projections and with , the following stability and approximation results are standard.

###### Lemma 3.1.

([27]) Let be an integer with . Then we have, for any and ,

 ∥v−Q0jv∥0,K+hK|v−Q0jv|1,K \apprlehsK|v|s,K,∀v∈Hs(K), ∥v−Q0jv∥0,∂K \apprlehs−1/2K|v|s,K,∀v∈Hs(K), ∥v−Qbjv∥0,∂K \apprlehs−1/2K|v|s,K,∀v∈Hs(K), ∥Q0jv∥0,K ≤∥v∥0,K,∀v∈L2(K), ∥Qbjv∥0,e ≤∥v∥0,e,∀v∈L2(e).

By using the trace theorem, the inverse inequality, and scaling arguments metioned in [27], we can get the following lemma.

###### Lemma 3.2.

For all , , and , we have

 ∥w∥0,~q,∂K\apprleh−1~qK∥w∥0,~q,K+h1−1~qK|w|1,~q,K,

In particular, for all ,

 ∥w∥0,~q,∂K\apprleh−1~qK∥w∥0,~q,K.
###### Lemma 3.3.

([6]) Let . For all and , the following estimates hold:

 ∥∇vh0∥0,K \apprle∥∇w,mvh∥0,K+h−1/2K∥Qblvh0−vhb∥0,∂K, (3.1) ∥∇w,mvh∥0,K \apprle∥∇vh0∥0,K+h−1/2K∥Qblvh0−vhb∥0,∂K, (3.2)

We introduce the following semi-norms: for any ,

 \interleavevh\interleave2:= ∥∇w,mvh∥20+∥τ1/2(Qblvh0−vhb)∥20,∂Tfh, ∥qh∥2:= ∥qh0∥20+∑K∈Tfh∥∇w,kqh∥20,K, \interleavesh\interleave2:= ∥∇w,msh∥20+∥τ1/2(Qblsh0−shb)∥20,∂Th.

Here we recall that . It is easy to see that the above three semi-norms are norms on , and , respectively (cf. [6]). In addition, from the lemma above it follows

 ∥∇hvh0∥0\apprle\interleavevh\interleave,∀vh∈V0h. (3.3)
###### Remark 3.1.

We note that the estimates (3.1), (3.2), and (3.3) also hold for all due to the fact that .

###### Lemma 3.4.

([14]) For all , there exists an interpolation such that

 ∑K∈Th∥sh0−Iksh0∥20,K \apprle∑e∈εhhe∥[sh0]∥20,e, ∑K∈Th∥∇(sh0−Iksh0)∥20,K \apprle∑e∈εhh−1e∥[sh0]∥20,e.

From this lemma it follows that, for all , there exists an interpolation such that

 ∑K∈Th∥vh0−Ikvh0∥20,K \apprle∑e∈εhhe∥[vh0]∥20,e, (3.4) ∑K∈Th∥∇(vh0−Ikvh0)∥20,K \apprle∑e∈εhh−1e∥[vh0]∥20,e. (3.5)
###### Lemma 3.5.

For all and , we have

 ∥vh0∥0,~q ≤C~q1\interleavevh\interleave, (3.6) ∥sh0∥0,~q ≤C~q2\interleavesh\interleave, (3.7)

where when , when , and , are positive constants only depending on .

###### Proof.

For all , we apply the Sobolev embedding theorem and Poincáre inequality to get

 ∥Ikvh0∥0,~q\apprle∥Ikvh0∥1\apprle∥∇Ikvh0∥0. (3.8)

From (3.5), (3.3), the definition of , and the projection property of , it follows

 ∥∇Ikvh0∥0 \apprle∥∇hvh0∥0+(∑e∈εhh−1e∥[vh0]∥20,e)12 (3.9) \apprle\interleavevh\interleave.

Using the Sobolev embedding theorem and the inverse inequality once again, by the properties of the projection-mean operator ([27]) and the fact that when and when , we have

 ∥vh0−Ikvh0∥0,~q ≤∥vh0−Πhvh0∥0,~q+∥Πhvh0−Ikvh0∥0,~q \apprleh∥∇hvh0∥0,~q+∥Πhvh0−Ikvh0∥1,2 \apprleh1−(d2−d~q)∥∇hvh0∥0,2+∥vh0−Πhvh0∥1,2+∥vh0−Ikvh0∥1,2 \apprle∥∇hvh0∥0+∥∇h(vh0−Ikvh0)∥0 \apprle\interleavevh\interleave,

which, together with (3.8) and (3.9), yields the desired estimate (3.6).

Similarly, we can obtain (3.7). This finishes the proof. ∎

For any nonnegative integer and any , we introduce the local Raviart-Thomas(RT) element space

 RTj(K)=[Pj(K)]d+xPj(K).

Lemmas 3.6-3.8 show some properties of the projection which can be founded in ([5].Page 9-10).

###### Lemma 3.6.

For any , implies .

###### Lemma 3.7.

For any and , there exists a unique such that

 ⟨PRTjv⋅ne,wj⟩e =⟨v⋅ne,wj⟩e,∀wj∈Pj(e),e∈∂K, (3.10) (PRTj