A high-order discontinuous Galerkin method for natural convection problems

TUDelft,TheNetherlands,2006

A HIGH-ORDER DISCONTINUOUS GALERKIN METHOD

FOR NATURAL CONVECTION PROBLEMS

F. Bassi, A. Crivellini

Dipartimento di Ingegneria Industriale

Universitàdi Bergamo

viale Mar oni 5,24044 Dalmine (BG),Italy

e-mail: fran es o.bassiunibg. it, andrea. rivelliniunibg. it

Keywords: In ompressibleow,Navier-Stokesequations,Dis ontinuousGalerkinmethod,

Arti ial ompressibility,Riemann solver, Boussinesq approximation

Abstra t. Dis ontinuous Galerkin (DG) methods have proved to be very well suited for

the onstru tionof robust high-order numeri al s hemesonunstru tured andpossibly non

onforming grids for a wide variety of problems. In this paper we onsider natural

on-ve tion ow problems and present a high-order DG method for their numeri al solution.

The governing equations are the in ompressible Navier-Stokes (INS) with the Boussinesq

approximationto represent buoyan y ee ts and the energy equationto des ribe the

tem-perature eld.

The method here presented is an extension to natural onve tion ows of a novel

high-order DG method for the numeri al solution of the INS equations, re ently proposed in

Bassi et al. 1

. The distinguishing feature of this method is the formulation of the invis id

interfa e ux whi h is based on the solution of lo al Riemann problems asso iated with

the arti ial ompressibility perturbation of the in ompressible Euler equations. The

dis- retization of the vis ous term follows the well established DG s heme named BR2 2,3

.

Themethodisfully impli itand thesolution isadvan edin timeusing eithera rstorder

ba kward Euler or a se ond order Runge-Kutta s heme.

Toassessthe apabilitiesoftheDGmethodpresentedinthispaperwe omputedse ond-,

third- andfourth-order spa e-a uratesolutions of several ben hmark problems onnatural

onve tion in two-dimension al avities.

1 INTRODUCTION

Dis ontinuous Galerkin(DG)methodshaveprovedtobesuitedforthe onstru tionof

robusthigh-ordernumeri als hemesonarbitraryunstru turedandpossiblynon

onform-ing grids for a wide variety of problems. The appli ation of the DG spa e dis retization

to in ompressible ows has been onsidered in a few re ent works. Liu and Shu 4

in a series of papers 57

Co kburn and oworkers proposed and thoroughly analyzed the

Lo al Dis ontinuous Galerkin (LDG) method for the Stokes, Oseen and Navier-Stokes

equations.

One of the key ingredientsof DG methodsis the formulation of interfa e (numeri al)

uxes, whi hprovide weak ouplingof dis ontinuous solutions inneighbouring elements.

In the invis id ompressible ase the numeri al ux is often omputed by exploiting the

hyperboli nature of the equations as the (approximate or exa t) solution of a Riemann

problem. Inthe in ompressible ase, however, theequationsare nolongerhyperboli due

tothela kofthe unsteadyterm inthe ontinuityequationand itisthereforenotpossible

to ompute the interfa e ux followingthe same approa h.

In this paper we ontinue the development of the DG method for the INS equations

introdu ed in Bassi et al. 1

The novelty of this method is the formulation of the invis id

numeri aluxwhi hisbasedonthesolutionoftheRiemannproblemforthe

in ompress-ibleEulerequationswitharelaxedin ompressibility onstraint. Itisimportanttoremark

that the idea of relaxing the in ompressibility onstraint by an arti ial ompressibility

term is employed only forthe onstru tionof the interfa e uxes.

The numeri al experiments reported in Bassi et al. 1

show that the onvergen e rate,

using polynomials of degree

k

for all the variables, is

k + 1

for the velo ity omponents and at least

k

for the pressure. In that paper it is also shown that the method is well suited for INS omputationsathigh Reynolds number.

Here the method is extended to solve the INS equations oupled with the energy

onservation equation using the Boussinesq approximation inthe body for e term of the

momentum equation. This set ofequations isasuitablephysi almodel todes ribemany

natural onve tion problems. In thefollowingse tionswewillshowthat thetreatmentof

the additionalenergy equation inthe Riemannproblemis straightforward.

As regards the DG spa e dis retization of the diusive terms we have used the well

established DG s heme introdu ed inBassi et al. 2

Finally,thespa edis retizedequationsareintegratedintimebymeansoffullyimpli it

rst andse ond ordera urate methodsand theresultinglinearsystems are solved using

LU dire t orGMRES iterative solvers.

The DGmethodheredeveloped has beenvalidatedby omputingtwosteady 810

(with

linear and non linear buoyan y ee ts) and one unsteady 8,9,11

thermally driven avity

ow problems.

2 GOVERNING EQUATIONS

We onsider the onservation equations ofmass, momentum and energy foran

in om-pressible uid. Followingthe Boussinesq approximation, buoyan y ee ts are a ounted

for by introdu ing a variable density body for e term in the momentum equation. The

density inthis term depends onthe temperature a ording tothe equation

where

ρ

is the density,

ρ

r

is the referen e density,

∆T

is the dieren e between the temperatureand its referen e value and

β

is the thermal expansion oe ient.

The non dimensional form of the governing equations is obtained assuming a length

s ale

L

r

, a temperature dieren e s ale

∆T

r

and the buoyan y velo ity s ale

U

r

=

√

gL

r

β∆T

r

, where

g

isthe a elerationof gravity. With the above hoi es the governingequations read:

∇

_{· u = 0,}

∂u

∂t

+ ∇ · F(u, p) − ∇ ·



Pr

√

Ra · Pr

∇

_u



= θj,

∂θ

∂t

+ ∇ · g(u, θ) − ∇ ·



1

√

Ra · Pr

∇

_θ



= 0,

(2) in

[0, t] × Ω

,

Ω ⊂ R

N

, where

N ∈ {2, 3}

is the number of spa e dimensions,

j

is the unit ve torinthe

y

dire tion,

u

and

θ

denotethedimensionlessvelo ityve torandtemperature dieren e,

Pr

is the Prandtl numberand

Ra = Pr

gL

3

r

β∆T

r

ν

2

,

the Rayleighnumber(

ν

is the kinemati vis osity), isthe main parameter hara terizing the ow eld.

The onve tive terms of the governing equations are writtenin onservation form

be- ausethis isthe naturalstartingpointforthe invis iduxdis retizationdes ribed inthe

followingse tion. The uxes

F(u, p)

and

g(u, θ)

are given by

F(u, p)

def

_{= u ⊗ u + pI = u}

i

u

j

+ pδ

ij

,

g(u, θ)

def

= uθ = u

i

θ,

where

p

is the dimensionless deviation from hydrostati pressure and

I

is the identity matrix.

3 DISCONTINUOUS GALERKIN DISCRETIZATION

In this se tionthe DG dis retizationof the governing equationsis rst introdu edfor

the invis id part of Eq. (2), and then ompleted with the dis retization of the vis ous

terms. The invis idandvis ous formsof Eq.(2)are looselydenoted asEuler and

e

x

y

P

n

−

n

+

K

+

K

−

Figure1: Normalsandlo alframeat quadraturepoint

P

onedge

e

.

3.1 Euler equations

Negle ting the diusive terms,the weak formof Eq. (2) reads:

−

Z

Ω

∇

_{q · u dx +}

Z

∂Ω

qu · n dσ = 0,

Z

Ω

v

_·

∂u

∂t

dx −

Z

Ω

∇

_v

_{: F(u, p) dx +}

Z

∂Ω

v

_{⊗ n : F(u, p) dσ =}

Z

Ω

v

_{· jθ dx,}

Z

Ω

q

∂θ

∂t

dx −

Z

Ω

∇

_{q · g(u, θ) dx +}

Z

∂Ω

qg(u, θ) · n dσ = 0,

(3)

for arbitrarytest fun tions

v

and

q

.

In order to onstru t a DG dis retizationof Eq. (3),we onsider atriangulation

T

h

=

{K}

of an approximation

Ω

h

of

Ω

, that is we partition

Ω

h

intoa set of non-overlapping elements

K

(not ne essarily simpli es). We denote with

E

0

h

the set of internal element fa es, with

E

∂

h

the set of boundary elementfa es and with

E

h

= E

0

h

∪ E

h

∂

their union. We moreover set

Γ

0

h

=

[

e∈E

0

h

e,

Γ

∂

_h

=

[

e∈E

∂

h

e,

Γ

h

= Γ

0

h

∪ Γ

∂

h

.

(4)

Thesolutionisapproximatedon

T

h

asapie ewisepolynomialfun tionpossibly dis ontin-uousonelementinterfa es,i.e.weassumethefollowingspa esettingsfortheapproximate

solution

u

h

,

p

h

and

θ

h

:

u

h

∈ V

h

def

= [V

h

]

N

,

p

h

, θ

h

∈ Q

h

def

= V

h

,

(5) with

V

h

def

=



v

h

∈ L

2

(Ω) : v

h

|

_K

∈ P

k

(K) ∀K ∈ T

h

for some polynomialorder

k ≥ 1

,being

P

k

(K)

the spa e of polynomials of global degree atmost

k

onthe element

K

.

Inordertosimplifythepresentation,itis onvenienttointrodu esometra eoperators

whi h generalize those dened in Arnold et al 3

. On a generi internal fa e

e ∈ E

Figure1, we dene

[[v]]

def

= v

+

⊗ n

+

+ v

−

_{⊗ n}

−

,

[[q]]

def

= q

+

n

+

+ q

−

n

−

,

(6)

where

v

denotes a generi ve tor quantity and

q

a generi s alar quantity. Noti e that

[[v]]

is a tensor quantity, and

[[q]]

is a ve tor quantity, i.e. this operator always in reases the tensor rankby one.

Moreover, we introdu e the average operator

{·}

def

=

(·)

+

+ (·)

−

2

,

(7)

whi happliestos alars,ve torsortensors. These denitions anbesuitablyextended to

fa es interse ting

∂Ω

a ounting for the weak impositionboundary onditions as shown inŸ 3.3.

The dis rete ounterpart of Eq. (3) for a generi element

K ∈ T

h

then reads:

−

Z

K

∇

_h

_q

_h

_{· u}

_h

_{dx +}

Z

∂K

q

h

u

h

|

_K

· n dσ = 0,

Z

K

v

h

·

∂u

h

∂t

dx −

Z

K

∇

_h

_v

_h

_{: F(u}

_h

_{, p}

_h

_{) dx}

+

Z

∂K

v

_h

_{⊗ n : F(u}

_h

_|

_K

, p

h

|

K

) dσ =

Z

K

v

_h

_{· jθ}

_h

dx,

Z

K

q

h

∂θ

h

∂t

dx −

Z

K

∇

_h

_q

_h

_{· g(u}

_h

_{, θ}

_h

_{) dx +}

Z

∂K

q

h

g(u

h

|

K

, θ

h

|

K

) · n dσ = 0.

(8)

To introdu e a ouplingbetween the degrees of freedom belonging to adja ent elements

andtoensure onservation,wesubstitutetheuxes

u

h

|

K

,

F(u

h

|

K

, p

h

|

K

)

and

g(u

h

|

K

, θ

h

|

K

)

withsuitablydenednumeri aluxes

u(u

b

±

h

, p

±

h

)

,

F(u

b

±

h

, p

±

h

)

and

b

g(u

±

h

, p

±

h

, θ

±

h

)

. Weremark that the stability and a ura y properties of the method strongly depend on the hoi e

of su h numeri aluxes.

SummingEq.(8)overtheelementsweobtaintheDGformulationofproblem(3)whi h

then requires tond

u

h

∈ V

h

and

p

h

, θ

h

∈ Q

h

su hthat

for all

v

h

∈ V

h

and

q

h

∈ Q

h

.

The key idea adopted to ompute

b

u

,

F

b

and

b

g

is toredu e the problemof ux ompu-tation tothe solutionof a planarRiemannproblemas inthe ompressible ase. Inorder

to re over the hyperboli hara ter of the equations, the in ompressibility onstraint is

relaxed by adding an arti ial ompressibility term to the ontinuity equation. At ea h

quadrature point

P

on

Γ

h

we thereforesolvethe Riemannproblemfor the equations

1

c

2

∂p

∂t

+

∂u

∂x

= 0,

∂u

∂t

+

∂ (u

2

_{+ p)}

∂x

= 0,

∂v

∂t

+

∂uv

∂x

= 0,

∂θ

∂t

+

∂uθ

∂x

= 0,

(10)

with initial datum

(u, p, θ) =

(

(u

−

_h

, p

−

_h

, θ

−

_h

)

if

x < 0,

(u

+

_h

, p

+

_h

, θ

+

_h

)

if

x > 0,

where

x

denotes a lo ally dened axis oriented as the normal ve tor

n

+

pointing out of

K

+

and lo ated insu h a way that

x = 0

at

P

,see Figure1.

Denoting with

(u

∗

, p

∗

, θ

∗

)

the solution of the Riemannproblemon the spa e-time line

x/t = 0

, we nallyset

b

u(u

±

h

, p

±

h

) = u

∗

,

F(u

b

±

_h

, p

±

_h

) = F(u

∗

, p

∗

),

b

g(u

±

_h

, p

±

_h

, θ

±

_h

) = g(u

∗

, θ

∗

).

It is worth noti ingthat, unlikein Co kburn et al. 57

, we donot split the numeri al ux

related to the adve tion term from that related to the pressure term, thus establishing

a stronger link between pressure and velo ity whi h seems to enhan e the properties of

the method. Due to spa e limitations, the pro edure for the determination of the state

(u

∗

, p

∗

, θ

∗

)

willnotbereportedhere. WerefertoBassietal. 1

forathoroughdes riptionof

thesolutionoftheRiemannproblemfor

(u

∗

, p

∗

)

. Asregardsthesolutionoftheadditional energy equation in system (10), we simply observe that the equation for

θ

has exa tly the same form as that for the tangential velo ity omponent

v

. Hen e, the form of the solutionfor

v

∗

reported in Bassi etal.

1

holdsalso for

θ

∗

. 3.2 Navier-Stokes equations

Many te hniques are available for the DG spa e dis retization of the diusive terms:

a omplete survey an be found in Arnold et al. 3

In this paper we hoose the form, rst

proposed in Bassi et al. 2

governingequations: nd

u

h

∈ V

h

and

p

h

, θ

h

∈ Q

h

su h that

−

Z

Ω

h

∇

_h

_q

_h

_{· u}

_h

_{dx +}

Z

Γ

h

[[q

h

]] · b

u

(u

±

h

, p

±

h

) dσ = 0,

Z

Ω

h

v

h

·

∂u

h

∂t

dx −

Z

Ω

h

∇

_h

_v

_h

_{: (F}

_v

_(∇

_h

_u

_h

_{, r) + F (u}

_h

_{, p}

_h

_{)) dx}

+

Z

Γ

h

[[v

h

]] :



b

F

_v

(∇

h

u

±

h

, r

±

e

) + b

F

(u

±

h

, p

±

h

)



dσ =

Z

Ω

h

v

_h

_{· jθ}

_h

dx,

Z

Ω

h

q

h

∂θ

h

∂t

dx −

Z

Ω

h

∇

_h

_q

_h

_{· (g}

_v

_(∇

_h

_θ

_h

_{, r) + g (u}

_h

_{, θ}

_h

_{)) dx}

+

Z

Γ

h

[[q

h

]] · bg

v

(∇

h

θ

h

±

, r

±

e

) + bg(u

±

h

, p

±

h

, θ

±

h

)

dσ = 0,

(11)

for all

v

h

∈ V

h

and

q

h

∈ Q

h

. The diusiveuxes inthe aboveequations are dened as:

F

v

(∇

h

u

h

, r)

def

= −

√

Pr

Ra · Pr

(∇

h

u

h

+ r ([[u

h

]])) ,

b

F

v

(∇

h

u

±

_h

, r

±

e

)

def

= −

√

Pr

Ra · Pr

({∇

h

u

h

} + η

e

{r

e

([[u

h

]])}) ,

g

_v

(∇

h

θ

h

, r)

def

= −

√

1

Ra · Pr

(∇

h

θ

h

+ r ([[θ

h

]])) ,

b

g

v

(∇

h

θ

h

±

, r

±

e

)

def

= −

√

1

Ra · Pr

({∇

h

θ

h

} + η

e

{r

e

([[θ

h

]])}) ,

(12)

where

r

e

is the lifting operator, that an be dened for the jumps of both ve tor and s alar quantities, and

r

is given by

r

_(·)

def

=

X

e∈E

h

r

e

(·) ,

(13)

being

e

a generi fa e. When applied to the jump of the velo ity ve tor

u

,

r

e

is dened as the solutionof the followingproblem:

Z

Ω

h

r

e

(ψ) : τ

h

dx = −

Z

e

ψ

_{: {τ}

h

} dσ, ∀τ

h

∈ [V

h

]

N

2

, ψ ∈



L

1

(e)



N

2

,

(14)

whereas applied tothe jumpof the s alar quantity

θ

itis dened as:

Z

Ω

h

r

e

(w) · v

h

dx = −

Z

e

w

_{· {v}

h

} dσ, ∀v

h

∈ [V

h

]

N

, w ∈



L

1

(e)



N

.

(15)

It is possible to nd lowerbounds for the parameter

η

e

∈ R

+

ensuring stability of the

method. We referthe interested reader toArnold et al. 3

3.3 Boundary onditions

The DG dis retization is best suited for a weak enfor ement of boundary onditions.

This an easily be a hieved by properly dening a boundary state whi h, together with

the internal state,allows to ompute the numeri aluxes and the liftingoperator onthe

portion

Γ

∂

h

of the boundary

Γ

h

.

The wall-type boundary onditions have been implemented by dening the boundary

stateontheexteriorofboundaryfa esas

u

b

h

= −u

+

h

,

∇

u

b

h

= ∇u

+

h

and

p

b

h

= p

+

h

. Diri hlet boundary onditions for

θ

pres ribe

θ

b

h

and set

∇

θ

b

h

= ∇θ

+

h

, while for adiabati wallswe set

θ

b

h

= θ

+

h

and

∇

θ

b

h

= −∇θ

+

h

. For all these ases the external boundary state exa tly repla es

u

−

h

,

θ

−

h

and

p

−

h

in the jumpoperators, in the numeri al uxes and in the lifting operator.

3.4 Time dis retization and linear system solution

All integralsappearing in Eq. (11) are omputed by means of Gauss integrationrules

with a number of integration points suited for the required a ura y. Cheaper

non-produ t formulae are preferred to tensor produ t ones when available. The quadrature

formulaearetakenfromtheen y lopædiaof ubatureformulaedevelopedandmaintained

by Cools 12

.

The dis rete problem orresponding toEq. (11) an bewritten as:

M

dW

dt

+ R (W) = 0,

(16)

where

W

is the global ve tor of unknown degrees of freedom and

M

is the global blo k diagonalmass matrix. La kingthetime derivativeofpressure inthegoverningequations,

the blo ks of

M

orresponding to the pressure degrees of freedom are identi ally zero. Eq. (16) denes a system of (nonlinear) ODEs whi h has been dis retized by means of

the se ond-order impli itRunge-Kutta s heme proposed by Iannelli and Baker 13 :

W

n+1

_{− W}

n

= Y

1

K

1

+ Y

2

K

2

,



M

∆t

+ α

∂R (W

n

₎

∂W



K

1

= −R (W

n

) ,



M

∆t

+ α

∂R (W

n

₎

∂W



K

2

= −R (W

n

+ βK

1

) ,

(17) where

∂R (W

n

_)/∂W

is the Ja obianmatrix of the DG spa e dis retizationand the

on-stants

α

,

β

,

Y

1

,

Y

2

are given by

α =

2 −

√

2

2

,

β = 8α



1

2

− α



,

Y

1

= 1 −

1

8α

,

Y

2

= 1 − Y

1

.

Noti e thatthe one-step ba kwardEuler andCrank-Ni olson s hemes an beobtained

Eq. (17) requires to solve a linear system of the form

Ax

+ b = 0

. However, sin e the matrix

A

is the same for the two steps, the Ja obian matrix needs tobe evaluated only on e.

The matrix

A

an be regarded as an

N

K

× N

K

blo k sparse matrix where

N

K

is the number of elements in

T

h

and the rank of ea h blo k is

N

K

DOF

× (N + 2)

, being

N

K

DOF

the numberof degrees of freedom for ea h variable in the generi element

K

. Thanks to the DG dis retization here adopted the degrees of freedom of a generi element

K

are only oupled with those of the neighbouring elements and the number of nonzero blo ks

for ea h (blo k) row

K

of the matrix

A

is therefore equal to the number of elements surrounding the element

K

plus one.

The Ja obianmatrix ofthe DGdis retizationhas been omputedanalyti allywithout

any approximation and, using very large time steps, the method an therefore a hieve

quadrati onvergen e in the omputation of steady state solutions. For the ba kward

Euler s hemeand inthe limit

∆t → ∞

Eq. (17)is infa tidenti al toone iterationof the Newton method appliedto the steady dis rete problem.

TosolveEq.(17)our ode an useone ofthenumerous approa hes (dire toriterative,

sequential or parallel) available in the PETS 14

library (Portable Extensible Toolkit for

S ienti Computations),the software upon whi h our DG ode relies for the purpose of

parallelization. The parallelization is based on grid partitioning a omplished by means

of the METIS 15

pa kage. In this ontext ea h pro essor owns the data related to its

lo al portion of the grid and the data on remote pro essors are a essed through MPI,

the standard formessage-passing ommuni ation. Thanks tothe ompa tness ofour DG

method onlythe data owned by the nearneighborelementsatpartitionboundaries need

tobeshared among dierent pro esses.

4 NUMERICAL RESULTS

Inthisse tionwe onsider threeben hmark problemsofnatural onve tionowwithin

en losed re tangular avities. The boundary onditionspres ribe thetemperatureonthe

verti al walls and no heat ux on top and bottom walls. In all ases the gravity for e is

dire ted downwards.

4.1 First ben hmark problem

The rst ben hmark problem isthe well known steady thermally driven avity ow.

The ow domainis asquare

(0, 1)

2

lled witha Boussinesq uid (air with

Pr = 0.71

). The problem, with Rayleigh numbers up to

10

6

, was sele ted as the test ase of a

numeri alben hmark in1983 16,17

. More re ently several resear hers extended the

ben h-mark in luding the

Ra = 10

7

and

Ra = 10

8

ases 9,1820

. Being very lose to the riti al

Rayleighnumber beyond whi hthe owbe omesunsteady (

(1.82 ± 0.01) · 10

8

,a ording

toLe Quéréetal.), the

Ra = 10

8

ase represents asevere test forany numeri almethod.

Here we present the numeri al results for the

Ra = 10

7

and

Ra = 10

8

Time steps

R

es

idu

al

s

10

20

30

40

50

60

10

-15

10

-13

10

-11

10

-9

10

-7

10

-5

10

-3

10

-1

10

1

(a) LUdire tsolver,

40

× 40

grid

Time steps

R

es

idu

al

s

20

40

60

80

100

120

140

160

10

-15

10

-13

10

-11

10

-9

10

-7

10

-5

10

-3

10

-1

10

1

(b) GMRES(60)iterativesolver,

80

× 80

grid

Figure2: Firstben hmarkproblem(

Ra = 10

8

),typi al onvergen ehistory.

ea hRayleighnumberwe have omputed

P

1

,

P

2

and

P

3

solutionsontwo stru turedgrids with

40 × 40

and

80 × 80

quad elements, lustered near the walls. The maximum grid spa ing,

h

max

= 0.045

, is the same for the two grids, whilst the minimum grid spa ing is

h

min

= 0.012

for the oarse grid and

h

min

= 0.0012

for the ne grid. The results here reported have been omputed with the arti ial ompressibility parameter

c = 1

.

The steady solutions of this problem have been omputed using the rst order

ba k-wardEuler s hemewithprogressivelyhighertimestepsasthesteadystateisapproa hed.

Even for the highest Rayleigh number solutions onverged to steady state in few

itera-tions. Typi al onvergen e histories of residuals onboth grids using dire t and iterative

linearsystemsolversare displayedinFigure2. Usually,higher-ordersolutionsare started

from the lower-order ones and this pra ti e turns out very useful to redu e the overall

omputational ost of higher-order omputations.

Figure3 shows the temperatureand the velo ity omponents isolines and the

qualita-tiveimprovement of the omputed solutionsin reasing the degree of polynomial

approx-imation.

A quantitative omparison of our most a urate solution (

80 × 80

grid,

P

3

elements) withthose reportedby otherauthors an befound inTables1 and2. The resultsused to

omparethesolutionsaretheNusseltnumber(average,maximumandminimumandtheir

lo ations) along the hot wall and the velo ity omponents (maxima and their lo ations)

along the enterlines. Following Le Quéré and Gjesdal maxima and minima have been

found by evaluating the polynomial approximation at

1001

equidistant points and the velo ity omponentshavebeen s aledwith

U = Pr

√

(a)

θ

,

P

1

elements (b)

θ

,

P

3

elements

( )

u

,

P

1

elements (d)

u

,

P

3

elements

(e)

v

,

P

1

elements (f)

v

,

P

3

elements

Figure 3: Firstben hmarkproblem(

Ra = 10

8

Present Mayne 20 Gjesdal 9 LeQuéré 18 work

Nu

16.523 16.387 16.523 16.523

Nu

max

39.394 41.025 38.940 39.395

y

0.018 0.039 0.020 0.018

Nu

min

1.366 1.380 1.366 1.366

y

1.0 1.0 1.0 1.0

u

max

· 10

2

4.699 4.594 4.699 4.699

y

0.879 0.884 0.879 0.879

v

max

· 10

2.211 2.222 2.211 2.211

x

0.021 0.021 0.021 0.021 Table 1: First ben hmark problem (

Ra = 10

7

),

80

_{× 80}

grid,

P

3

elements. Present Mayne 20 Gjesdal 9 LeQuéré 18 work

Nu

30.225 29.626 30.229 30.225

Nu

max

87.224 91.209 86.700 87.236

y

0.008 0.007 0.009 0.008

Nu

min

1.919 2.044 1.919 1.919

y

1.0 1.0 1.0 1.0

u

max

· 10

2

3.219 2.831 3.219 3.219

y

0.928 0.946 0.928 0.928

v

max

· 10

2.222 2.223 2.222 2.222

x

0.012 0.013 0.012 0.012 Table 2: First ben hmark problem (

Ra = 10

8

),

80

_{× 80}

grid,

P

3

elements.

been evaluated as:

Nu = (∇

h

θ

h

+ r ([[θ

h

]])) |

w

· n,

wherethesubs ript

w

meansthatthe fun tioninsideparenthesesis omputedatthe wall and the gradient of the temperature in ludes the ontribution due to interfa e

dis onti-nuities. The data reported in Tables 1 and 2 indi ate that our results are in ex ellent

agreementwith the best solutions available inthe literature.

Table 3 summarizes the results of a onvergen e study for the

Ra = 10

8

ase. In

this Table we have alsoreportedinside parentheses the data omputedwith the arti ial

ompressibility parameter

c = 0.1

. On the whole the results onrmthe a ura y of the methodhere proposed. The inuen e ofthe arti ial ompressibility parameteris almost

negligible for the high order approximation while larger hanges an be observed for

P

1

solutions. Thisis trueinparti ularforthe maximumof

u

velo ity omponentas an also be seen in Figure 4(a). Numeri al experiments with

c = 10

resulted in a poor or even impossible onvergen e of the solution.

4.2 Se ond Ben hmark problem

This ben hmark problem, proposed by Mi halek et al. 10

,is aslightly modied version

of the previous test ase. The square avity is in this ase lled with water and the

temperatures of the isothermal walls, i.e.

T

h

= 10

◦

C and

T

c

= 0

◦

C, are around the

temperature of maximum density. The anomalous thermal variation of water density is

des ribed by the polynomialfun tion

ρ (T ) = 999.840281167108 + 6.73268037314653 · 10

−2

T

− 8.94484552601798 · 10

−3

_T

2

+ 8.78462866500416 · 10

−5

_T

3

− 6.62139792627547 · 10

−7

_T

4

,

(18)

40

× 40

40

× 40

40

× 40

80

× 80

80

× 80

80

× 80

LeQuéré 18

P

1

P

2

P

3

P

1

P

2

P

3

Nu

29.971 30.177 30.221 30.234 30.225 30.225 30.225 (30.028) (30.196) (30.224) (30.238) (30.225) (30.225)

Nu

max

81.707 90.401 88.961 86.915 87.231 87.224 87.236 (87.260) (92.259) (88.875) (86.880) (87.226) (87.224)

y

0.011 0.011 0.011 0.009 0.008 0.008 0.008 (0.000) (0.011) (0.010) (0.009) (0.008) (0.008)

Nu

min

1.865 1.904 1.928 1.934 1.919 1.919 1.919 (1.749) (1.925) (1.927) (1.920) (1.919) (1.919)

y

1.0 0.998 1.0 1.0 1.0 1.0 1.0 (1.0) (0.998) (1.0) (1.0) (1.0) (1.0)

u

max

· 10

2

2.50646 3.14109 3.21221 2.72554 3.20706 3.21888 3.21875 (2.64660) (3.18459) (3.21558) (3.02003) (3.21979) (3.21929)

y

0.930 0.931 0.928 0.933 0.927 0.928 0.928 (0.914) (0.931) (0.928) (0.926) (0.927) (0.928)

v

max

· 10

2.27924 2.23883 2.22349 2.21639 2.22226 2.22238 2.22239 (2.29722) (2.24012) (2.22361) (2.21657) (2.22224) (2.22238)

x

0.012 0.012 0.012 0.012 0.012 0.012 0.012 (0.000) (0.012) (0.012) (0.012) (0.012) (0.012)

Table3: First ben hmark problem (

Ra = 10

8

), omparison of solutions with dierent polynomial

ap-proximationandgridresolution.

The Prandtlnumberisset to

13.31

andthe Rayleighnumber, atthereferen e temper-ature

T

r

= 0

◦

C, is equal to

1.503 · 10

6

. It isknown that the variation of vis osity is not

negligible in the range of temperatures here onsidered, but for simpli ity and with the

aimof buildinganeasily omparable ben hmark solution allthe uid hara teristi s are

assumed to be onstant.

The omputations have been performed on equally spa ed, su essively rened (from

4 × 4

to

128 × 128

)quad grids, using the same time integrationmethodsof the previous problem.

As suggested in Mi halek 10

, omputed solutions have been ompared on the basis of

the average Nusselt number on the hot wall and of the maxima and minima of velo ity

omponents along enterlines. Values and positions of maxima and minima have been

evaluated from the solution interpolated at 1001 equally spa ed points. The results of a

grid onvergen e studyusing

P

3

elementsare summarized inTable 4,wherethevelo ities havebeens aledwith

U = Pr/ (νL

r

)

. For omparisonpurposesthelast olumninTable4 reports the referen e solution omputed by Mi halek et al.

10

on a

301 × 301

grid using the nite dieren e ode FRECONV3V.

y

u

0

0.2

0.4

0.6

0.8

1

-0.04

-0.03

-0.02

-0.01

0

0.01

0.02

0.03

0.04

P

1

, c =0.1

P

1

, c =1.0

P

3

, c =0.1

P

3

, c =1.0

(a)

u

along

x = 0.5

, inuen e of the arti ial ompressibilityparameter,

80

× 80

grid

x

v

0

0.02

0.04

0.06

0.08

0.1

-0.05

0

0.05

0.1

0.15

0.2

0.25

P

1

P

2

P

3

(b)

v

along

y = 0.5

, dierent polynomial ap-proximations,

40

× 40

grid

Figure4: Firstben hmarkproblem(

Ra = 10

8

),velo ity omponentsalong enterlines.

with two ounter-rotating vorti es in the zone where the temperature rea hes the

maxi-mum density value. One an alsoappre iatethat thisowpatternisalready aptured on

the

8 × 8

oarse grid. Finally,the velo ity proles of Figure6 show that the

P

3

solution is almostgrid onverged onthe

16 × 16

grid.

4.3 Third ben hmark problem

Unlike the previous test ases the last ben hmark problem is unsteady. In this ase

the uid is air (

Pr = 0.71

), the omputationaldomain is an en losed 8:1 (height:width) re tangular avity and the Rayleighnumber, assumingthe widthas referen e length

L

r

, is

3.4 · 10

5

.

This problemwasproposed attheMIT Conferen e onComputationalFluid andSolid

Dynami s 8

in2001. Several ontributors provided 32 sets of solutiondata using abroad

range of omputational methods. The results obtained by Le Quéré at al. 11

using a

Chebyshev ollo ation spe tral method on

48 × 180

points were sele ted as the baseline for the omparison and the ranking of the methods. In a re ent paper Gjesdal et al.

9

reported the results of a spe tral element simulation(

18

th

order ona

4 × 20

grid) whi h fully agree with LeQuéré's data.

In this work the problemhas been omputed using

P

1

,

P

2

and

P

3

polynomial approxi-mationsand the grid resolutions suggested atthe MITConferen e, namely

21 × 101

and

(a)

θ

,

8

× 8

grid (b)

θ

,

128

× 128

grid

( )

u

,

8

× 8

grid (d)

u

,

128

× 128

grid

(e)

v

,

8

× 8

grid (f)

v

,

128

× 128

grid

y

u

0

0.2

0.4

0.6

0.8

1

-100

-50

0

50

100

4x4

8x8

16x16

32x32

64x64

128x128

(a)

u

along

x = 0.5

x

v

0

0.2

0.4

0.6

0.8

1

-200

-100

0

100

200

300

4x4

8x8

16x16

32x32

64x64

128x128

(b)

v

along

y = 0.5

Figure6: Se ondben hmarkproblem,velo ity omponentsalong enterlines,

P

3

elements.

4

× 4

8

× 8

16

× 16

32

× 32

64

× 64

128

× 128

T.Mi halek 10

P

3

P

3

P

3

P

3

P

3

P

3

Nu

6.0073 6.4492 6.4660 6.4679 6.4680 6.4680 6.47

u

min

· 10

2

-1.4209 -1.3420 -1.3123 -1.3098 -1.3098 -1.3098 -1.31

x

0.751 0.692 0.708 0.708 0.709 0.709 0.71

u

max

19.7656 4.3480 2.8749 2.8893 2.8903 2.8903 2.94

x

0.375 0.354 0.376 0.382 0.382 0.382 0.38

v

min

· 10

2

-1.3054 -1.7201 -1.7524 -1.7515 -1.7514 -1.7514 -1.75

x

0.705 0.691 0.705 0.705 0.705 0.705 0.70

v

max

· 10

2

1.7962 2.2130 2.1744 2.1712 2.1709 2.1709 2.17

x

0.064 0.047 0.041 0.042 0.042 0.042 0.04

u

min

· 10

-7.7215 -8.0671 -7.7815 -7.78276 -7.7827 -7.7827 -7.78

y

0.304 0.305 0.286 0.286 0.286 0.286 0.29

u

max

· 10

8.5708 8.5147 8.6746 8.6788 8.6789 8.6789 8.67

y

0.885 0.892 0.893 0.892 0.892 0.892 0.89

v

min

· 10

-11.0830 -10.133 -8.7888 -8.7872 -8.7856 -8.7856 -8.77

y

0.251 0.273 0.261 0.261 0.261 0.261 0.26

v

max

16.4307 6.3343 6.4743 6.4918 6.4919 6.4920 6.48

y

0.630 0.626 0.647 0.649 0.648 0.649 0.65

ratio equal to 1.2 and 1.05, respe tively.

Sin e the problemis unsteady the solutionhas been advan ed in time using the more

a urate se ond-order Runge-Kutta time stepping s heme. In our fully impli it time

dis retization,the time stepsize is notsubje ttostabilitylimitationsand has been xed

onthebasis ofa ura yrequirements. From numeri alexperimentswefound thatatime

step equal to

τ /40

, where

τ = 3.4115L

r

/U

r

is the period of the solution omputed by LeQuéré,was su iently smalltoassure a urate results. Infa t, themaximum relative

dieren e between two

P

2

solutions on the

21 × 101

grid, omputed doubling the time step size from

τ /40

to

τ /20

, was less than

1%

.

Also for this unsteady problem the

P

k

solutionwas started fromthe

P

k−1

approxima-tion. Wefoundthatthe

P

1

approximationwasnotadequateto aptureowunsteadiness, even onthe nestgrid. Startingfromthe steady

P

1

results,the

P

2

and

P

3

solutionshave been advan ed in time sequentially. For ea h approximationwe found that about 15,000

time steps were largely su ient to rea h almostperfe tlyperiodi solutions.

21

× 101

21

× 101

21

× 101

41

× 201

41

× 201

LeQuéré 11 Gjesdal 9

P

2

P

2

P

3

P

2

P

3

∆t = 1/20τ

∆t = 1/40τ

∆t = 1/40τ

∆t = 1/40τ

∆t = 1/40τ

θ

1

· 10

1

2

.6546

2

.6550

2

.6548

2

.6548

2

.6547

2

.6548

2

.6561

θ

′

1

· 10

2

4

.2443

4

.2046

4

.2872

4

.2350

4

.2885

4

.2689

4

.2774

u

1

· 10

2

5.6177

5.6175

5.6343

5.6262

5.63421

5.6345

5.6356

u

′

1

· 10

2

5

.4505

5

.3903

5

.5035

5

.4323

5

.5049

5

.4767

5

.4880

p

14

· 10

3 (

1

)

−1.8618

−1.8508

−1.8492

−1.8654

−1.8515

−1.8536

-p

′

14

· 10

2 (

1

)

2.0279

2.0087

2.0443

2.0226

2.0448

2.0355

-ω

1

−2.4030

−2.4012

−2.3695

−2.3637

−2.3721

-

−2.3678

ω

′

1

1

.0802

1

.0643

1

.0802

1

.0755

1

.0859

-

1

.0801

Nu

−4.57939

−4.57940

−4.57946

−4.57944

−4.57946

−4.57946

−4.57933

Nu

′

_{· 10}

3

_7.0750

_7.0075

_7.1241

_7.0504

_7.1241

_7.0918

_7.1026

b

u · 10

1

₂

_.3949

₂

_.3948

₂

_.3951

₂

_.3950

₂

_.3951

-

2

.395

b

u

′

_{· 10}

5

_3.4021

_3.3366

_3.3742

_3.3444

_3.3770

-

3.354

b

ω

3.0166

3.0165

3.0171

3.0177

3.0171

-

3.0171

b

ω

′

_{· 10}

3

₃

_.1964

₃

_.1612

₃

_.2078

₃

_.1773

₃

_.2090

-

3

.198

Table5: Third ben hmarkproblem, onvergen estudy.

Table 5 summarizes our omputational results. All mean data reported in the Table

havebeenaveragedover

20

periods. The monitored quantitiesarethe averagedvalueand the peak to peak amplitude of the variables at the rst ontrol point (pres ribed at the

1

p

14

is the pressure dieren e between the rst and the fourth ontrol point (

x

4

= 0

.819Lr

and

MIT onferen e as

x

1

= 0.181L

r

and

y

1

= 7.37L

r

), the averaged Nusselt number along the hot wall and the averaged velo ity magnitude (

u

b

) and vorti ity (

ω

b

) over the whole omputational domain. Our results are in ex ellent agreement with the referen e data

andwementionthat,usingtheoverallmetri softhe errorsdened attheMIT onferen e

for the omparison of the methods, both the oarse and the ne grid

P

3

solutions would have been ranked rst, i.e. losest to the referen e solution, and the

P

2

solutions sixth and se ond.

All the solutions of this test ase have been omputed running the ode in parallel

and using the blo k Ja obi pre onditioned GMRES algorithm of the PETS library for

the linear algebra. A relative residual error of

10

−8

was used as stopping riterion for

the iterativesolver. Typi ally, the

P

3

omputationson the oarse grid took 120 GMRES iterationsto onverge the solution within the pres ribed toleran e at ea htime step.

Figure7: Thirdben hmarkproblem,isolinesoftheinstantaneoustemperaturedeviationsfrom thelo al

5 CONCLUSIONS

Inthisworkanovelhigh-orderDGmethodforthenumeri alsolutionofthe

in ompress-ibleNavier-Stokesequationshasbeen extendedtodealwithnatural onve tionproblems.

The idea of using the solution of perturbed 1D Riemann problems for the purpose of

inter-element ux omputation has been su essfully applied to the dis retization of the

onve tiveterm in the energy onservation equation.

The method has been validated by omputing a number of well known steady and

unsteady natural onve tion ben hmark problems. In all ases our DG solutionswere in

ex ellent agreement with the best referen e results for su h problems.

REFERENCES

[1℄ F.Bassi,A.Crivellini,D.A.DiPietro,andS.Rebay.Anarti ial ompressibilityux

forthedis ontinuousGalerkinsolutionofthein ompressibleNavier-Stokesequations.

J. Comput. Phys. ,arti le ele troni ally published, 2006.

[2℄ F. Bassi,S. Rebay, G.Mariotti,S. Pedinotti, and M. Savini. A high-order a urate

dis ontinuousniteelementmethodforinvis idandvis ousturboma hineryows.In

R. De uypere and G. Dibelius, editors,Pro eedings of the 2nd European Conferen e

onTurboma hineryFluidDynami sandThermodynami s,pages99108,Antwerpen,

Belgium,Mar h 57 1997. Te hnologis hInstituut.

[3℄ D.N. Arnold,F.Brezzi, B.Co kburn,and D.Marini. Uniedanalysisof

dis ontinu-ousGalerkinmethodsforellipti problems. SIAMJ.Numer.Anal.,39(5):17491779,

2002.

[4℄ J.-G. Liu and C.-W. Shu. A high-order dis ontinuous Galerkin method for 2D

in- ompressible ows. J. Comput. Phys. ,160:577596, 2000.

[5℄ B. Co kburn, G. Kans hat, D. S hötzau, and C. S hwab. Lo al dis ontinuous

Galerkin methods for the Stokes system. SIAM J. Numer. Anal., 40(1):319343,

2002.

[6℄ B. Co kburn, G. Kans hat, and D. S hötzau. The lo al dis ontinuous Galerkin

method for the Oseen equations. Math. Comp. , 73(246):569593, 2003.

[7℄ B. Co kburn, G. Kans hat, and D. S hötzau. A lo ally onservative LDG method

for the in ompressibleNavier-Stokesequations. Math. Comp. , 74:10671095, 2005.

[8℄ M. A. Christon, P. M. Gresho, and S. B. Sutton. Computational predi tability of

time dependent natural onve tion (in luding a ben hmark solution). Internat. J.

[9℄ T. Gjesdal, C. E. Wasberg, and B. A. P. Reif. Spe tral element ben hmark

simula-tionsof natural onve tionintwo-dimensional avities. Internat.J.Numer. Methods

Fluids, inpress, 2005.

[10℄ T. Mi halek, T. A. Kowalewski, and B.

ˇ

Sarler. Natural onve tion for anomalous

density variation of water: numeri al ben hmark. Progress in Computational Fluid

Dynami s ,5:158170, 2005.

[11℄ S. Xin and P. Le Quéré. An extended Chebyshev pseudo-spe tral ben hmark for

the 8:1 dierentiallyheated avity. Internat. J.Numer. Methods Fluids, 40:981998,

2002.

[12℄ R.Cools. AnEn y lopædiaof CubatureFormulas. J.Complexity ,19:445453,2003.

[13℄ G. S. Iannelli and A. J. Baker. A stiy-stable impli it Runge-Kutta algorithm for

CFD appli ations. AIAA Paper 88-0416,AIAA,1988.

[14℄ Satish Balay, Kris Bus helman, William D. Gropp, Dinesh Kaushik, Matthew G.

Knepley, Lois Curfman M Innes, Barry F. Smith, and Hong Zhang. PETS Web

page, 2001. http://www.m s.anl.gov/pets .

[15℄ G.KarypisandV.Kumar. METIS,asoftwarepa kage forpartitioningunstru tured

graphs,partitioningmeshes,and omputingll-redu ingorderingsofsparsematri es.

Te hni al Report Version 4.0, University of Minnesota, Department of Computer

S ien e/Army HPCResear hCenter, 1998.

[16℄ G.DeVahlDavis. Natural onve tionofairinasquare avity: aben hmark

numer-i alsolution. Internat. J. Numer.Methods Fluids, 3:249264, 1983.

[17℄ G. De Vahl Davis and I. P. Jones. Natural onve tion of air in a square avity: a

omparison exer ise. Internat. J. Numer. Methods Fluids,3:227248, 1983.

[18℄ P.LeQuéré.A uratesolutionstothesquarethermallydriven avityathighRayleigh

number. Comput.& Fluids, 20(1):2941,1991.

[19℄ D. C. Wan, B. S. V. Pantnaik, and G. W. Wei. A new ben hmark quality solution

for the buoyan y-driven avity by dis rete singular onvolution. Numeri al Heat

Transfer. Part B, 40:199228,2001.

[20℄ D. A. Mayne, A. S. Usmani, and M. Crapper. h-adaptive nite element solution of

highRayleighnumberthermallydriven avityproblem. Internat.J.Numer.Methods