A pressure-based solver for low-Mach number flow using a discontinuous Galerkin method

Delft University of Technology

A pressure-based solver for low-Mach number flow using a discontinuous Galerkin

method

Hennink, Aldo; Tiberga, Marco; Lathouwers, Danny

DOI

10.1016/j.jcp.2020.109877

Publication date

2021

Document Version

Final published version

Published in

Journal of Computational Physics

Citation (APA)

Hennink, A., Tiberga, M., & Lathouwers, D. (2021). A pressure-based solver for low-Mach number flow

using a discontinuous Galerkin method. Journal of Computational Physics, 425, [109877].

https://doi.org/10.1016/j.jcp.2020.109877

Important note

To cite this publication, please use the final published version (if applicable).

Please check the document version above.

Copyright

Other than for strictly personal use, it is not permitted to download, forward or distribute the text or part of it, without the consent of the author(s) and/or copyright holder(s), unless the work is under an open content license such as Creative Commons. Takedown policy

Please contact us and provide details if you believe this document breaches copyrights. We will remove access to the work immediately and investigate your claim.

This work is downloaded from Delft University of Technology.

Green Open Access added to TU Delft Institutional Repository

'You share, we take care!' - Taverne project

https://www.openaccess.nl/en/you-share-we-take-care

Otherwise as indicated in the copyright section: the publisher

is the copyright holder of this work and the author uses the

Dutch legislation to make this work public.

Contents lists available atScienceDirect

Journal

of

Computational

Physics

www.elsevier.com/locate/jcp

A

pressure-based

solver

for

low-Mach

number

flow

using

a

discontinuous

Galerkin

method

Aldo Hennink

∗

,

Marco Tiberga,

Danny Lathouwers

DepartmentofRadiationScienceandTechnology,DelftUniversityofTechnology,TheNetherlands

a

r

t

i

c

l

e

i

n

f

o

a

b

s

t

r

a

c

t

Articlehistory:

Received8November2019

Receivedinrevisedform22September 2020

Accepted26September2020 Availableonline2October2020 Keywords:

Low-Mach Variableproperties DiscontinuousGalerkin Pressurecorrection

Over thepasttwodecades,therehasbeenmuchdevelopment indiscontinuousGalerkin methodsforincompressibleflowsandforcompressibleflowswithapositiveMachnumber, but almost no attention has been paid to variable-density flows at low speeds. This paperpresentsapressure-baseddiscontinuousGalerkinmethodforflowinthelow-Mach number limit.Weuseavariable-density pressurecorrectionmethod,whichissimplified by solvingforthe massflux insteadofthevelocity.The fluidpropertiesdo not depend significantlyonthepressure,butmayvarystronglyinspaceandtimeasafunctionofthe temperature.

We payparticularattentiontothetemporaldiscretizationoftheenthalpyequation,and show that the specific enthalpy needs to be ‘offset’ with a constant in order for the temporalfinitedifferencemethodtobestable.Wealsoshow howonecansolve forthe specific enthalpyfrom the conservative enthalpy transport equation withoutneeding a predictorstepforthedensity.Thesefindingsdonotdependonthespatialdiscretization. Aseriesofmanufacturedsolutionswithvariablefluidpropertiesdemonstratefull second-ordertemporalaccuracy,withoutiteratingthetransportequationswithinatimestep.We alsosimulateaVonKármánvortexstreetinthewakeofaheatedcircular cylinder,and showgoodagreementbetweenournumericalresultsandexperimentaldata.

©2020ElsevierInc.Allrightsreserved.

1. Introduction

Severallow-speedflows ofpractical importance arecompressible, that is,the velocity isnot divergence-free.Thiscan occur duetomixing, ordueto atemperature-dependentdensitynearaheat source.Anexampleis heattransferin low-Machnumberflowsofsupercriticalfluids,whereallfluidpropertiesvarystronglywiththetemperature,butdonotdepend significantlyonthepressure.Mostflowsolversuseeitherapressure-basedapproachandassumeadivergence-freevelocity field, ora fullycompressible(density-based)formulation. Neither ofthesemethods isdirectlyapplicable tocompressible flowsinthelow-Machnumberlimit.

Density-basedsolvers canbeused tosimulatezero-Mach flowsby approximatingtheflowwithalow,non-zeroMach number(e.g.,[46],[11]).Thishasoftenbeenusedforheattransferinsupercriticalfluidsatlowspeeds(e.g.,[14],[4]).This isexpensiveforseveralreasons.First,thetemporaldiscretizationneedstoresolveacousticeffects,andtheresultinglinear systemstendtobeverystiff.Second,thesystemoftransportequationsissolvedinacoupledway,whichismoreexpensive than usingatime-splittingmethod,though theperformancemaybe improvedwithsuitable preconditioning[31].Finally,

*

Correspondingauthor.

E-mailaddress:a.hennink@protonmail.ch(A. Hennink). https://doi.org/10.1016/j.jcp.2020.109877

the fluid propertiesare evaluated asafunction oftwo thermodynamic variables (usuallythe densityandthevolumetric enthalpy),sothatasplineinterpolationcostsfarmorememory,thuscomplicatingmassivelyparallelcalculations[14].

ThereisalsosubstantialexperiencewithdiscontinuousGalerkin(DG)discretizationsforincompressibleflows.Theseare either based on theintroduction of artificialcompressibility (e.g.[12]),or they solve forthe pressure (e.g. [41,8,37,35]). The artificial compressibility methodcan be morethan second-order accurate in time, though it requires the systemof transport equationsto be solved in a coupled manner (e.g., [6,5]). By choosing entropy variables as the unknowns, the DGmethodcan alsobeformulated inageneralwayforbothcompressibleandincompressibleflows,atthecostofgreat complexity(e.g.,[34]).Thereis,however,almostnoliteratureonsolvingthelow-Machnumberequationswitha pressure-baseddiscontinuousGalerkinmethod.

The only previous work of which we are aware is by Klein et al. [24,25], who used a SIMPLE scheme to march the transportequationsforwardintime,iteratingtheequationswithineachtimestep.Thisrequiredunder-relaxationinorder for theiteration to converge. Theysolved for thevelocity, so that apredictor forthe densityis neededin thetemporal derivativeofthemomentumequation.

We avoid this by solving for the mass flux rather than the velocity. Another advantage of this approach is that the divergenceterminthecontinuityequationdoesnothavetobeweighedbythedensity,sothatthedivergencematrixdoes notdependonthedensity.Thismakesthetransportequationslesstightlycoupled,anditsimplifiesthepressurecorrection method,becausethepressurematrixisconstantforeachtimestep.

Thetemporaldiscretizationoftheenthalpyequationrequiresspecialcare.Theconvectivetermintheprimitiveenthalpy equation cannot beupwindedwithafinite elementmethodifthevelocityis notdivergence-free.Wethereforediscretize theenthalpyequationinconservativeform,althoughwesolveforthespecific enthalpy,whichisaprimitivevariable.This causes a complicationin thetemporal discretization,which features theunknown density ata newtime step.A similar issueoccursinmultispeciestransport,andNajmet al.[28] devisedawide-spreadtwo-stepiterativemethodtostabilizethe temporalscheme.Thishassubsequentlybeenadaptedtohandlethestrongpropertyvariationsinsupercriticalfluids[29].

Wepresentanalternativemethodthatdoesnotuseanypredictorsolvesoriterationstohandletheunknowndensityat anewtimestep.Weshowthat theerrorinourapproximationcanbemadenegligiblecomparedtotheerrorinthefinite differencescheme,andthatthemethodcanbemadestablebyoffsettingthespecificenthalpywithaconstant.

The paper is structured asfollows. Section 2 defines the mathematical problem. Section 3 and 4 explain the spatial discretizationandthetime-splittingscheme,payingparticularattentiontothepolynomialsolutionspaces(Section3.1)and thediscretizationoftheviscousstress(Section3.2.1),whichdifferslightlyfrompreviousliterature.Section5and6treatthe densitypredictorintheenthalpyequationingreatdetail,andillustratetheresultswithnumericalexamples;thesesections donot dependonthespatial discretization,andthey canberead independently.Section 7verifiesthenumericalmethod andtheimplementationforthecoupledtransportequationswithaseriesofprogressivelymorechallengingmanufactured solutions,andvalidatesourapproachbyreproducingexperimentaldataforflowpastaheatedcircularcylinder.Finally,we drawourconclusionsbasedontheresultsinSection8.

2. Governingequations

Thetransportequationsinthelow-Machnumberlimitare

∂

ρ

∂

t

+ ∂

kmk

=

0 , (1a)

∂

mi

∂

t

+ ∂

k

(

ukmi

)

= ∂

k

τ

ik

− ∂

ip

+

Fi, (1b)

∂

ρ

h

∂

t

+ ∂

k

(

mkh

)

= −∂

kqk

+

Q . (1c)

Heret isthetime,

ρ

isthedensity,u isthevelocity,m

:=

ρu is

themassflux,p isthepressure,h isthespecificenthalpy, andF and Q areknownexternalsources.AssumingaNewtonianfluid,thestresstensoris

τ

i j

=

μ



∂

iuj

+ ∂

jui

−

2 3

∂

kuk

δ

i j



. (2)

Theheatfluxis

qi

= −

k

∂

iT

= −

k cp

∂

ih , (3)

whereT isthetemperature.k isthethermalconductivity,andcpisthespecificheatcapacity.Thelastequalityistechnically

an approximation becauseit neglects thedependenceofthe temperatureon thepressure,butthisis highly accurate for low-Machnumberflows,even forstronglyvariable fluidpropertiesinsupercritical fluids[33].The fluidproperties(

ρ

,

μ

,

cp,andk)areafunctionofT ,butdonotdependsignificantlyonthepressure.

•

Dirichletboundaries,denotedby



D,onwhichthemassfluxandthetemperature(and,therefore,theenthalpyandthe fluidproperties)aregiven,thatis,mi

=

mDi andT

=

TD;

•

outflowboundaries,denotedby



N,where

(

τ

ik

−

p

δ

ik

)

nk

=

fiNandknl

∂

lT

=

qNareprescribed.

3. Spatialdiscretization

This section details the spatial discretizationof the enthalpy andmass flux transportequations on a domain



with outwardnormaln andspacevectorr

= [

x

,

y

]

∈ 

.Thedomainismeshedintoasetofelements

T

.Wedenoteby

F

i,

F

D, and

F

N _{thesets ofinternalfaces, Dirichletboundaryfaces,andNeumannboundaryfaces. Thesetoffaces ofanelement}

T

∈

T

is

F

T,and

T

F isthesetofneighborsofface F .Eachinternal face F

∈

F

i hasanormalvectornF withanarbitrary

butfixeddirection.Thejumpandaverageoperatorsaredefinedas



x



:=

x+

−

x− and

{

x

} :=

1 2



x−

+

x+



, (4) where x±

(

r

∈

F

)

:=

lim ↓0x

(

r

∓



n F

₎

₍₅₎

indicate thevaluesatthetwosides. Onboundary elements,both thejumpandtheaverageequalthe internalnumerical value,andnF coincideswiththeoutwardnormalof



.

The solutionspacewithin each elementis spannedbyall polynomialsup toan order

P

. Wedenotethe orderofthe polynomialsfortheunknown X by

P

X.Forthenumericalexamplesinthispaper,thepolynomialorderisthesameonall

elements,thoughthisisnotarequirementofthenumericalmethod.

3.1. Solutionsspacesforenthalpyandpressure

In this paper the basis functions for the enthalpy and the pressure are always of the same order, that is,

P

h

=

P

p.

Neglectingthesourceterm,theenthalpyequationinitsconservativeformis

Dh

Dt

=

1

ρ

∂

l

(

k

∂

lT

)

, (6)

where D

/

Dt isthetotalderivative.If

∂

ρ

/∂

h

<

0 forallh,thenthisbecomes

1

ρ

D

ρ

Dt

=

∂

ρ

∂

h 1

ρ

2

∂

l

(

k

∂

lT

)

= −

β

ρ

cp

∂

l

(

k

∂

lT

)

= −

α

β

1 k

∂

l

(

k

∂

lT

)

, (7)

where

β

:= −(

1

/

ρ

)∂

ρ

/∂

T isthe thermalexpansibility,and

α

:=

k

/

ρc

p isthethermaldiffusivity. Comparethe continuity

equation

1

ρ

D

ρ

Dt

= −∂

kuk. (8)

Regardless ofthe Péclet number orthe compressibility, the right-hand-sides ofEqs. (7) and(8) canbe arbitrarily small inanypartof thedomain,inwhichcasethe enthalpyequation andthecontinuityequation are almostequivalent. Since thecontinuityequationinitsweakformisweighedby thepressurebasis functions,itseemsinconsistenttousedifferent polynomialordersforthepressureandtheenthalpyif

∂

ρ

/∂

h

<

0.

Furthermore,ournumericalexperiments(notshownhere)demonstratethat,forlowdiffusivity,thespatialdiscretization isunstableif

P

h

>

P

p,eveniftheenthalpy isapassivescalar,the densityisconstant, andwe simulatea laminarsteady

state.Thereasonfortheinstabilityisthattheconvectivevelocityintheenthalpydiscretizationmustsatisfythecontinuity equationuptotheorder

P

h,inaweaksense.ThisobservationcontrastswithKleinet al.[24],whousedapolynomialorder

forthetemperaturethatwashigherthanforthepressure.Thisdiscrepancyseems tobebecausetheirtestsweredone at alowPrandtlnumberof0.7.Wealsoobtainedaccurateresultswith

P

h

>

P

p forthe2Dbuoyancy-drivencavityatvarious

Richardson numbers,butthe instabilitymanifestsitself whenthethermaldiffusionvanishes.We willshow thisinfuture work.

3.2. Momentumdiscretization

Wesolveforthemassfluxm fromtheconservativetransportequation(1b).Kleinet al.[24] havepreviouslysolvedthe sameconservative equation bytakingthe velocity u asthe unknown,butthispresentsdifficulties. First,one wouldhave tohandlethestronglyfluctuatingdensityinthetemporalderivative

∂(

ρu

)/∂

t.Second,solvingforu wouldcomplicate the pressurecorrectionscheme,aswewillseeinSection4.Anotheralternativeistosolveforu fromtheprimitiveequation

ρ



∂

ui

∂

t

+

uk

∂

kui



= ∂

k

τ

ik

− ∂

ip

+

Fi, (9)

asincommonforincompressibleflows.(See[17,1] forexampleswithavariabledensity.)Unfortunatelywecannotupwind theconvectiveterminEq.(9),becauseitcannotberewrittenas

∂

k

(

uiuk

)

whenthecompressibility(

∂

kuk)isnon-negligible.

Ourapproachoftakingm astheunknown intheconservativeequation slightlycomplicatestheimplicittreatment ofthe viscousterm,whichislinearin

∂

iu,not

∂

im.Section3.2.1detailshowthiscanbehandledwithaDGmethod.

Thesemi-discreteweakformofthemomentumequationis

Find m

∈

V , such that, for all v

∈

V ,



 vk

∂

mk

∂

t

+

a conv

₍

_u

_,

_m

k

,

vk

)

+

avisc

(

m

,

v

)

=

lconv

(

u

,

mDk

,

vk

)

+

lvisc

(

v

)

−

adiv

(

v

,

p

)

+





Fkvk, (10)

where V isthe solution spaceof the massflux. The divergence operator adiv _{is (see, e.g., [41], [21,pp. 92], or [36, pp.} 250–252]) adiv

(

v

,

q

)

= −



T∈T



T q

∂

kvk

+



F∈Fi_∪FD



F

{

q

}



vk



nkF, (11)

whereq isascalarthatliesinthesolutionspaceofthepressure.

3.2.1. Discretizationoftheviscousstress

Toderiveadiscretizationfortheviscousterm,rewritetheviscousstressintermsofthemassfluxas

τ

=

Lvisc

₍

_m

₎

_,where

Lvisc_{i j}

(

m

)

=

μ

ρ



∂

imj

+ ∂

jmi

−

2 3

(∂

kmk

) δ

i j

−

dimj

−

djmi

+

2 3

(

dkmk

)δ

i j



, (12)

isalinearoperator,and

di

:=

1

ρ

∂

i

ρ

. (13)

Weuseageneralizationofthesymmetricinteriorpenalty(SIP)method,givenbythediscretebilinearoperator

avisc

(

w

,

v

)

=



T∈T



T Lvisc_kl

(

w

) ∂

lvk

+



F∈Fi_∪FD



F

η

F

_

_w k

 

vk



−



F∈Fi_∪FD



F





vk



L_klvisc

(

w

)

+



wk



Lvisc_kl

(

v

)



n_lF (14) and lvisc

(

v

)

=



F∈FD



F



η

FmD_k vk

−

mDk Lvisckl

(

v

)

nl



+



N f_kNvk. (15)

This reduces to the regular SIP method for constant-density, constant-viscosity flows when substituting

(

μ

/

ρ

)∂

lwk for

Lvisc_kl

(

w

)

.Comparedtoother interiorpenalty methodsandthelocalDGmethod,theadvantagesoftheSIPmethodarethe optimalconvergencerateforallpolynomialorders,itsadjointconsistency,anditscompactsmallstencil[18].

This discretization of the viscous term can be compared to what is usually done for compressible flows, where the systemofequations (1) issolved fora fullstate vector U

:= [

ρ

,

m

,

ρ

h

]

T. Inthat caseallelliptic termsinEqs. (1) canbe writtenas

∂

k

(

G

(

U

)∂

kU

))

,whereG

(

U

)

isahomogeneitytensorthatdoesnotcontainanygradientsoftheunknowns.(See,

e.g., [18,23].)This tensoris then kept fixed during an iteration step,while

∂

iU can be treatedimplicitly.Eq. (12) shows

that eachtermoftheviscousstresstensorisbilinear inU :they areproducts ofm and

ρ

.The lastthreetermscontaina gradientof

ρ

,andthereforethemassfluxisevaluatedinthehomogeneitytensor,andfrozenwithinaniterationstep.This differsfromourimplicitdiffusiondiscretization,where

ρ

isfrozenineachterminEq.(12),andm istreatedimplicitly.

Our approachalso differs fromthat ofKlein et al. [24], in that we treatall terms inthe viscous stress (Eq. (2)) ina time-implicitmanner, whereasthey onlydothisforthefirstterm(

μ

∂

iuj).Inourtreatment thevelocitycomponentsare

coupled.Wehavenoaprioriestimateforthedifferenceinmagnitudebetweentheeffectsofthefirstterm(

μ

∂

iuj)andits

transpose(

μ

∂

jui)ontheviscousforce

∂

i

τ

i j,especiallywhentheviscosityvariesstronglyinspace.Notethegradientsinthe

effectiveviscositywouldincreasegreatlyifalargeeddysimulation(LES)modelwereincluded. FollowingHillewaert[20,pp.30],wesetthepenaltyparameterto

η

F

₌

_max T∈TF



CTcard

(

F

T

)



F



_leb



T



leb



max T∈TF

(

K

|

T

)

, (16)

where K

=

μ

/

ρ

isthediffusionparameter,card

(F

T

)

isthenumberoffacesofelementT ,and

·leb

denotestheLebesgue

measure(whichisthelength,area,orvolumein1,2,or3dimensions).ThefactorCT dependsonthetypeofelementsin

themesh:forapolynomialorder

P

,CT

= (

P +

1

)

2 forquadrilateralsandhexahedra,CT

= (

P +

1

)(

P +

2

)/

2 ontriangles,

andCT

= (

P +

1

)(P +

3

)/

3 for tetrahedra.We compute thepenalty parameterpointwise, ratherthan averagingover the

face.

3.2.2. Discretizationofconvectiveterm

Thediscretizationoftheconvectivetermisgivenby

aconv

(

b

, φ, ψ )

= −



T∈T



T

ψ

bk

∂

k

φ

+



F∈Fi



F



ψ



HF

(

b

, φ)

+



N

(

nkbk

) ψ φ

+



D max



0

,

nkbDk



ψ φ

, (17) lconv

(

b

, φ, ψ )

= −



D min



0

,

nkbDk



ψ φ

. (18)

Here

φ

and

ψ

are scalars,b is theconvective field, and HF isthe numericalflux functionon aface F .Weusethe Lax-Friedrichsflux,givenby

HF

(

b

, φ)

=

1 2



φ



α

F

₍

_b

₎

₊

_nF k

{φ

bk

}

(19) and

α

F

₍

_b

₎

₌

_{f max} T∈TF





n_kFbk





T . (20)

Theconstant f is f

=

2 whenconvectingthemassflux,and f

=

1 whenconvectingascalar,suchastheenthalpy[9].Note thatweevaluate

α

F _{pointwiseintheintegralinEq.(17),ratherthanaveragingovertheface.}

ItiswellknownthatimposingaDirichletboundaryconditionforthevelocityisnotstableatanoutletofa convection-dominatedflow,sothatwewouldnormallyhavenkbD_k

<

0 on



D.HereweneverthelessincludetheDirichletboundaryterm

inaconv_{,becausewewilluseitintheTaylor-GreenvortexinSection7.2,asisstandardpracticeforthatlaminarbenchmark} case(e.g.,[13]).

3.3. Enthalpydiscretization

Wesolveforthespecificenthalpy h becauseitfixes thethermodynamicstate foragiventhermodynamicpressure pth. Thisis not trueforthevolumetricenthalpy (H

:=

ρh):

thepair

(

pth

,

H

)

doesnot uniquely determinethedensityor the otherfluidproperties.Choosingh astheunknownalsosimplifiesthetime-implicittreatmentoftheFourierheatflux.

Aswiththemassfluxequation,wesolvetheenthalpyequationinconservativeform(Eq.(1c)),becausewecannotdeal with the convective term(uk

∂

kh) inthe primitive transport equation when

∂

kuk

=

0.Given a solution space W for the

enthalpy,whichisalsothesolutionspaceforthepressure,thesemi-discreteweakformoftheenthalpyequationis

Find h

∈

W , such that, for all v

∈

W ,





v

∂(

ρ

h

)

∂

t

+

a

conv

₍

_m

_,

_h

_,

_v

₎

₊

_aSIP

₍

_h

_,

_v

₎

₌

_lconv

₍

_m

_,

_hD

_,

_v

₎

₊

_lSIP

₍

_v

₎

₊





Q v , (21)

whereaSIPandlSIParestandardSIPbilinearandlinearformstodiscretizetheFourierheatflux.TheSIPpenaltyparameter isasinEq.(16), withadiffusioncoefficient K

=

k

/

cp.Notethattheconvective discretizationisthesameasforthemass

flux,exceptthattheconvectivefieldism,ratherthanu

=

m

/

ρ

.(Thatis,b

=

m inEqs.(17)-(20).)Thisisanotheradvantage ofsolvingforthepair(m,h).

4. Pressurecorrectionmethod

Wediscretizethetemporalderivativesinboththeenthalpyandmomentumtransportequationswithstandard backward-difference formulae(BDF), therebyfollowing previousDG literature(e.g.,[41,24,26]).Forthemassflux thisis straightfor-ward:foraconstanttimestepsize

δ

t,

∂

m

∂

t

≈

γ

0

δ

tm n

₊

q



i=1

γ

i

δ

tm n−i_{, (22)}

Table 1

Coefficientsforthebackwarddifferenceformulaofvariousorders.

γ0 γ1 γ2 γ3

BDF1 1 −1

BDF2 3/2 −2 1/2

BDF3 11/6 −3 3/2 −1/3

where mn _{is the mass flux at time step n. The weights}

_{

_γ

i

}

qi=0 are listed inTable 1. The temporal discretization of the enthalpyequationrequiresmorecare,aswewillexplainlaterinSection5.

Apressure-correctionschemeisusedto splitthecontinuityandthemomentum equations,sothat theycanbe solved inasegregatedway.Thealgorithmtofindthesolutionvectorsmn_,hn_,pn _{atanewtimestepn isasfollows.}

1. Obtainpredictorsfor

ρ

∗,

(

k

/

cp

)

∗,andm∗ withasecond-orderextrapolationfromprevioustimesteps:

(

·)

∗

=

2

(

·)

n−1

− (·)

n−2. (23)

2. Solvefortheenthalpyhn atthenewtimestep,usingtheabovepredictorsforthediffusionconstant(

(

k

/

cp

)

∗)andthe

convectivefield(m∗).

LinearizetheimplicittimetermwiththepredictorsusingeitherofthemethodsthatareexplainedinSection5. 3. Findapredictorforthemassfluxfromalinearsystemthatcorrespondstothesemi-discreteweakform(10):



_γ

0

δ

tM

+

N



ˆ

m

= −

Dpn−1

+

f . (24)

Here M isthemassmatrix,N containstheimplicitpartsofthediffusionandconvectiondiscretizations,andf collects all explicitterms(i.e.,boundaryconditions,theexternalforce,andtheexplicittermsinthetemporalfinitedifference scheme).ThematrixN containstheviscosityandthedensity,whichareevaluatedatthenewtimestepasafunction ofhn.Theconvectivefieldisestimatedasm∗

/

ρ

n_.

4. SolvethepressurePoissonequation

δ

t

γ

0 A

δ

p

=

Dm

ˆ

−

G



∂

ρ

∂

t



n



, (25)

whereA isthepressurematrix,forwhichwederiveanexpressionbelow,and

G[·]

denotestheGalerkinprojectiononto thesolutionspace.Thetemporalderivativeofthedensityisestimatedwithasecond-orderfinitedifferencescheme:



∂

ρ

∂

t



n

≈

1

δ

t



3 2

ρ

n

₋

₂

_ρ

n−1

₊

1 2

ρ

n−2



. (26)

5. Correctthepressureandthemassflux:

mn

= ˆ

m

−

δ

t

γ

0

M−1DT

δ

p , (27)

pn

=

pn−1

+ δ

p . (28)

Thetest casesinSection 7.3willshowfull second-ordertemporalaccuracyintheenthalpy andthemassflux,evenifthe fluidpropertiesarenon-trivialfunctionsoftheenthalpy.

Tofindanexpressionforthepressurematrix A,useEq.(25) toeliminate

δ

p inEq.(27),andleft-multiplybyD,giving

Dmn

=

Dm

ˆ

−

D M−1DA−1



Dm

ˆ

−

G



∂

ρ

∂

t



n



. (29)

Clearlythissatisfiesthesemi-discretecontinuityequation

−

Dm

+

G

∂

ρ

∂

t



=

0 (30)

if A

=

ALDG

₌

_{D M}−1_D_{.As explained byShahbaziet al. [41], A}LDG _{iseffectivelya localDGdiscretization fora diffusion} operator,sowecanreplaceitbyanSIPdiffusionoperator

ALDG

≈

ASIP, (31)

Equal-order discretizations (i.e.,

P

p

=

P

m) are knownto be unstable forfinite element methods.Cockburn et al. [10]

addressedthisprobleminthecontextofaconstantkinematicviscosityandusingthe ALDG _{discretizationbyaddingaterm} oftheform



F∈Fi

γ

0



F



leb

ν



F



φ

 

q



(32)

forbasisandtestfunctions

φ

andq tothepressurematrix.Weset

γ

0

=

1,andadjusttheabovepenaltytermtoavariable viscositybytaking

ν

insidetheintegral,andtakingtheminimalvalueof

ν

onbothsidesoftheface.

Inthetestcasesthatfollow,wehavefoundnonoticeabledifferencebetweenusingLDGandSIPpressurediscretizations, which is in line withprevious findings (e.g., [25, pp. 33–45]). Interestingly, Shahbazi [40] has successfully used the SIP pressurematrixwithanequal-orderdiscretizationwithoutextrapressurestabilization[40,pp.48–65].Ourtests(notshown here)alsoindicatethat, forequal-orderdiscretizations withoutpressurestabilization,theLDGpressure matrixisunstable forallreasonabletimesteps,whereasusingtheSIPmatrixisfeasibleforawiderangeofpracticaltimestepsizes,though italwaysbecomesunstableinthelimit

δ

t

→

0.

Notethatsolvingforthemassflux simplifiesthepressure-correctionscheme.Ifwehadinsteadsolvedforthevelocity, thenthedensitywouldhadtohavebeenincorporatedintothedivergenceoperator D inEq.(25),andintothemassmatrix

M inEq.(24).

Thisextensionofthepressure-correctionmethodtocompressibleflowshassometimesprovedunstableinfinite differ-enceschemesthatwereappliedtomixingflowswithlargedensityratios(ofapproximatelymorethanafactorof3),because thecontinuityequationwasnotsatisfiedintheinviscidlimit;seeNicoud[30].Itisnotcertainwhetherthesameinstability wouldoccurforthediscontinuousGalerkinmethodpresentedhere;ourexperiencesofarhasnotexposedinstabilitieswith largedensityratios.Nicoudsuggestedadifferentgeneralizationofthepressurecorrectionmethodtovariable-densityflows, wherethedensityisincorporatedintothepressurematrix,ratherthanontheright-handsideofEq.(25).Thelarge advan-tage oftheapproachpresentedhereisthatthepressure matrixisthesameatalltime steps.Wecanthereforeassemble itonce,andprecomputetheincompleteCholeskypreconditionerforthelinearsolver.Furthermore,thecondition number ofthe diffusionmatrix A worsensifitincludes avariable coefficientthat dependsonthedensity.Forthesereasons,the pressuresolvesaremuchcheaperwithaconstantpressurematrix.

5. Linearizingthetemporalderivativeoftheenthalpy

Solving fora primitivevariable(h) withtheenthalpy equation inconservativeform slightlycomplicatesthetemporal derivative,becauseitisweighedbythetemperature-dependentdensity.Tostudythestabilityandconvergenceofthetime steppingscheme,weconsiderthespace-independententhalpyequation:

d

(

ρ

h

)

dt

= −λ

h

+

Q , (33)

where

λ

isaconstant,andQ

=

Q

(

t

)

.Usinganimplicitfinitedifferencescheme,theenthalpyandthecorrespondingdensity canbeestimatedatatimestepn by

γ

0

δ

t

(

ρ

h

)

n

₊

q



i=1

γ

i

δ

t

(

ρ

h

)

n−i

_{= −λ}

_hn

₊

_Qn_{. (34)}

Duetothevariabledensity,thisequationisnotlinearintheunknown hn_{.Wethereforeconsidertwolinearizationsinh}n_,

whichwetermmethod#1andmethod#2.

Both ofthese methodsuse apredictorh∗ anda corresponding

ρ

∗ that areclose tohn and

ρ

n_{.Thispredictor canbe}

obtained in several ways, such as by extrapolating from previous time steps. When solving the full system (1a)-(1c), a predictorfor

ρ

n _{canalsobeobtainedbysolving thecontinuityequation.Theanalysesinthissectionareforageneral(h}∗_,

ρ

∗),thoughwewillmakethereasonableassumptionthat

(

h∗

−

hn

)

isatleast

O(δ

t

)

. Thetwolinearizationmethodsareasfollows.

Method #1 isperhapsthemostobviousapproach:let

ρ

n

_≈

_ρ

∗_{,resultinginanapproximationh}[1]

_≈

_hn _{thatisgivenby}

γ

0

δ

t



ρ

∗h[1]



+

q



i=1

γ

i

δ

t

(

ρ

h

)

n−i

_{= −λ}

_h[1]

₊

_Qn_{. (35)}

Method #2 isbasedonaTaylorexpansionof

(

ρh

)

naboutthepredictor:

(

ρ

h

)

n

≈ (

ρ

h

)

∗

+



∂(

ρ

h

)

∂

h



_∗

_

hn

−

h∗



=



∂(

ρ

h

)

∂

h



_∗ hn

−



h2

∂

ρ

∂

h



_∗ . (36)

SubstitutingthisintoEq.(34) yieldsanapproximationh[2]

≈

hn,givenby

γ

0

δ

t



∂(

ρ

h

)

∂

h



_∗ h[2]

+

q



i=1

γ

i

δ

t

(

ρ

h

)

n−i

₌

γ

0

δ

t



h2

∂

ρ

∂

h



_∗

− λ

h[2]

+

Qn. (37)

Notethatmethod#1iseffectivelyasinglestepinafixed-pointiteration,whereasmethod#2isasinglestepinaNewton iteration.

5.1. Errorestimatesandstability

Theerrorsandthestability ofmethods#1and#2canbeanalyzedbyusingaTaylorseriesfor

ρ

n _{aboutthepredictor,}

thatis,

ρ

n

=

∞



k=0 1 k

!



∂

k

ρ

∂

hk



_∗



hn

−

h∗



k . (38)

Definethefollowingdeviationsfromthenon-linearfinitedifferenceequation(34):

error in the predictor:



∗

:=

h∗

−

hn, linearization error in method #1:



[1]

:=

h[1]

−

hn, linearization error in method #2:



[2]

:=

h[2]

−

hn.

(39)

Thederivationsaretedious,anddeferredtotheAppendix.Herewesummarizethemaintheoreticalresults.

Thefirstresultisanapriorierrorestimate.AppendixAshowsthatbothEq.(35) andEq.(37) canberewrittenas

γ

0

δ

t

(

ρ

h

)

n

₊

q



i=1

γ

i

δ

t

(

ρ

h

)

n−i

_{= −λ}

effhn

+

Qeffn , (40) where

λ

eff

= λ +

O





∗2

/δ

t



and Q_effn

=

Qn

+

O





∗2

/δ

t



for method #1, (41a)

and

λ

eff

= λ +

O





∗3

/δ

t



and Q_effn

=

Qn

+

O





∗3

/δ

t



for method #2. (41b)

Thatis,theapproximationsinmethod#1and#2areequivalenttotheoriginalEq.(34),exceptthat

λ

andQn arereplaced bytheireffectivevalues,whicharerelatedtotheerrorinthepredictor.

Asecondimportantresultregardsthestabilityofthelinearizationmethods.InAppendixA.1weshowthattheerrorfor method#1isrelatedtotheerrorinthepredictoras



ρ

∗

+

δ

t

γ

0

λ





[1]

₌



−



∂

ρ

∂

h



_∗ hn





∗

+

O





∗2



. (42)

Notethat



[1] _{vanishesuptofirstorderin}

_

∗_ashn

_→

_{0.Eq.(42) alsoshowsthatmethod#1cannotalwaysbemadestable}

byiteratingwithinatimestep,thatis,bycalculatinganewpredictor

ρ

∗ fromtheestimateh[1],andrepeatingEq.(35).For thestabilityoftheiteration,theerrorinthenewapproximationneedstobesmallerthantheerrorinthepredictor,thatis,

|



[1]

_|

_<

_|

_

∗

_|

_{.Takingthelimit}

_δ

_t

_→

_{0 inEq.(42),thisconditionisonlymetif}



hn



<





ρ

∂

ρ

/∂

h





∗

=



cp

β



_∗ , (43)

wherewehavemadethereasonableassumptionsthat



∗ isatleast

O (δ

t

)

,andthat

∂

ρ

/∂

h

<

0.NotethatEq.(43) cannot becheckedduringacalculation,becauseitdependsontheunknownhn.

Similarly,iteratingmethod#2withinatimestepisstableif

|



[2]

_|

_<

_|

_

∗

_|

_{.AppendixA.2showsthat}



∂(

ρ

h

)

∂

h



_∗

+

δ

t

γ

0

λ





[2]

=



∂

ρ

∂

h



_∗

+

1 2



∂

2

ρ

∂

h2



∗ hn





∗2

+

O





∗3



. (44)

Sincewe canreasonablyexpectthat theerrorinthepredictor(



∗) isatleast

O(δ

t

)

,wealwayshave

|



[2]

_|

_<

_|

_

∗

_|

_(i.e.,the



∂(

ρ

h

)

∂

h



∗

=

0 . (45)

Inotherwords,thevolumetricenthalpy

(

ρh

)

mustbeastrictlymonotonicfunctionofthespecificenthalpyh.

This restriction formethod #2also follows moredirectly fromEq. (37) in the limit ofsmall time steps, becausethe coefficientofh[2] _{cannotvanish.Inpracticeonewillwanttosatisfythestrongerrelation}

∂(

ρ

h

)

∂

h

>

0 (46)

toensurethattheenthalpydiscretizationispositivedefinite.

Weconjecturethatthestabilityrequirements(Eq.(43) formethod#1;Eq.(45) formethod#2)mustalwaysbesatisfied inthelimitofsmalltimesteps,evenwhenthelinearizationisnotiteratedwithinatimestep.Itseemsreasonabletoexpect thatastablenumericalmethodcanbeiteratedwithoutdiverging.ThisissupportedbythenumericaltestsinSection6

5.2. Properscalingoftheenthalpyequation

Curiously,theanalysesintheprevioussubsectionhaveledtostabilityrequirements(Eqs.(43),(45))thatdependonthe fluid properties, andthey are not satisfied forall fluids. For example,the volumetricenthalpy in supercritical fluidscan eitherincreaseordecreasewiththetemperatureduetothestrongthermalexpansion,therebyviolatingEq.(45).

Thisproblemcanbeaddressedbysolvingforadifferentvariable

˜

h

:=

h

−

h0. (47)

Eq.(1c) thenbecomes

h0R

+

∂(

ρ

h

˜

)

∂

t

+ ∂

l



mlh

˜



= ∂

l



k cp

∂

lh

˜



+

Q , (48)

where R

:= ∂

ρ

/∂

t

+ ∂

kmk

=

0 istheresidualofthecontinuityequation(1a).Thus h satisfies

˜

thesametransportequation

as h, andit can thereforebe discretized in the same way.This does not affect the convective or diffusive terms (since

∂

ih

˜

= ∂

ih),butitdoesmakeadifferenceforthetemporalderivative.

Inparticular,thestabilityrequirement(Eq.(43))formethod#1becomes



˜

hn



 ₌



hn

−

h0



<



cp

β



∗ . (49)

There isnoapriori valueforh0 that guaranteesstability,though itseems thath0 isbestchosen suchthat h

≈

h0 atthe averagetemperature.

Conversely, we can find an a priori lower bound for h0 when using method#2. The stability requirement(Eq. (46)) becomes

∂(

ρ

h

˜

)

∂ ˜

h

=

∂(

ρ

(

h

−

h0

))

∂

h

=

ρ

+ (

h

−

h0

)

1 cp

∂

ρ

∂

T

>

0

⇔

h0

>

h

−

cp

β

, (50)

and so method #2 can be made unconditionally stable and SPD by choosing h0 sufficiently large. In particular, if the temperatureisknowntolieinarange

[

Tmin

,

Tmax

]

,thenthetheoreticalminimumvalueis

hmin₀

=

max Tmin_≤T≤Tmax



h

−

cp

β



, (51)

which can of course be negative. Fig. 1 showsan example of the rescaled volumetric enthalpy

(

ρ

(

h

−

h0

))

forvarious choicesofh0forarealfluid.Thereisapracticallimitonthemagnitudeofh0,becausewearesolvingatransportequation for

ρ

(

h

−

h0

)

,whichbecomesequivalenttothedensity

ρ

forverylargevaluesof

|

h0|.

Ofcoursearescalingoftheunknowns,suchasinEq.(47),iscommoninCFDliterature,butitisusuallypresentedasa mere numericalconvenience (e.g.,[33]).Theaboveanalysesshow thattheaccuracyandstabilityofthenumericalscheme depend criticallyon aproper choiceofh0.Inpractice thismayrequiresome trialanderror,though theseanalyses offer usefulguidelines.

Fig. 1. Rescaledvolumetricenthalpy(ρ(h−h0))ofcarbondioxideatthesupercriticalpressureof8MPa,asafunctionofthetemperatureforvarious choicesofh0.Thefunctionismonotonicforh0≥hmin0 .Thethermodynamicreferencepointisplacedat(1bar,0◦C).Thedataarebasedon[45],accessed throughtheCoolPropsoftwarelibrary[7].(Forinterpretationofthecolorsinthefigure(s),thereaderisreferredtothewebversionofthisarticle.)

6. Testcaseforthespace-independentEnthalpyEquation1

Before solving the full systemoftransport equations,we clarify thetheoretical results forthe space-independent en-thalpy equation in Section 5with a numericalexample that isbased on a manufactured solution.Omitting the unitsof measurement,theexacttemperatureis

Tex

(

t

)

=

0

.

5

+

0

.

1 sin

(

2

π

t

)

(52)

with0

≤

t

≤

1.Theequationofstateis

ρ

=

ρ

0T

+

ρ

1

(

1

−

T

)

, (53)

andthespecificheatcapacityiskeptconstant,sothat

h

=

cpT

−

h0. (54)

Therequiredsourceterm Q

(

t

)

followsfromEq.(33).Forthenumericaltestswe let

ρ

0

=

0

.

5,

ρ

1

=

2,cp

=

1,and

λ

=

0

.

1.

The results presentedherewere all obtained at

δ

t

=

2−11 _{to investigatethe limit ofsmall time steps.We have checked} that loweringthetime stepsizeto

δ

t

=

2−14 doesnotaffectwhetherthenumericalmethodisstable. Tomakesurethat rounding errors donot play a significant role withthesetiny time steps, all calculationsinthis section were performed with128-bitfloatingpointprecision.

ThenumericalschemesweretestedwithvariousordersoftheBDFtimesteppingscheme.Thepredictorh∗ isobtained withan sth_{-orderextrapolationfromprevioustime steps(denoted byEXs),andthecorresponding}

_ρ

∗ _{isdeterminedfrom} theequationofstate.Specifically,

h∗

=

s



i=1

α

ihn−i

=

hn

+

O



δ

ts



. (55)

TheweightsareinTable2.

There aretwo numericalerrorsineach timestep: (i)theBDF error,whichisinherentinthe finitedifferencescheme, and(ii)thelinearizationerroringoingfromEq.(34) toeitherEq.(35) orEq.(37) whenusingmethod#1ormethod#2.If theEXs coefficientsareusedtoobtainapredictor,thentheerrorinthepredictoris



∗

:=

h∗

−

hn

=

O



δ

ts



.Formethod#1, Eq.(41a) thenimpliesalinearizationerrorof

O





∗2

/δ

t



=

O



δ

t2s−1



_{.ABDFq schememakesan}

_O



_δ

_tq+1



_errorpertime step,sothattheoverallorderofaccuracyismin

(

2s

−

1

,

q

+

1

)

.Similarly,Eq.(41b) impliesthattheoverallerrorpertime stepformethod#2isofordermin

(

3s

−

1

,

q

+

1

)

.Theorderofextrapolationshouldthereforesatisfy

s

≥



(

q

+

2

)/

2 for method #1,

(

q

+

2

)/

3 for method #2, (56)

1 _{Wehavepublishedthecodeforthefinitedifferencemethodforthespace-independententhalpyequationonGitHub[19].Itcanbeusedtoreproduce} theresultsinthissection.

Table 2

Coefficientsfor extrapolationfromprevioustimesteps. (SeeEq.(55)). α1 α2 α3 α4 EX1 1 EX2 2 −1 EX3 3 −3 1 EX4 4 −6 4 −1 Table 3

Orderofextrapolationfortheenthalpypredictor(Eq.(55))forthe linearizationmethodsdescribedinSection5.Theminimumvalues satisfyEq.(56);fromthemaximumvalueonward,Eq.(56) holds withstrictinequality.

Finitedifference coefficients

Method #1 Method #2 min max min max

BDF1 EX2 EX2 EX1 EX2

BDF2 EX2 EX3 EX2 EX2

BDF3 EX3 EX3 EX2 EX2

or else the linearizationerror dominates,and theusual orderof convergenceof the BDFscheme cannot be achieved. If strictinequalityinEq.(56) issatisfied,thenthelinearizationerrorisnegligible,andincreasingtheorderofextrapolationis pointless.Table3liststherangeofreasonableextrapolationorders.

Fig. 2showstheerror inthenumerical temperatureasafunction ofh0, byusing method#1 (Eq.(35)).The stability requirementinEq.(43) cannotbeverifiedduring thecalculation,becauseitdependsontheunknownhn (i.e.,thesolution to the finite difference scheme without linearizationerror; see Eq.(34)). But hn _{can be approximated by the numerical}

solutionh[1]_{,sothatthestabilityrequirementisapproximately}





h[1]



 _<



cp

β



_∗

.

(57)

AccordingtoTable3,thelinearizationerrorisnegligibleforextrapolationordersofatleast2,3,and3fortheBDF1,BDF2, andBDF3schemes.Forthesecasestheerror

|

h[1]

−

hn

|

isnegligible,andthenumericalschemeconvergesifandonlyifEq. (57) issatisfiedatalltimesteps.Notethattherangeofstablevaluesforh0 decreaseswithhigher-orderextrapolations,but allsimulationsconvergewithh0

=

0

.

5,whichisthevalueforwhichh isclosesttozero.

Fig.3showstheequivalenterrorplotsformethod#2(Eq.(37)).Forthecurrentequationofstate,wehaveanexplicit,a prioriexpressionforthestabilitycriterioninEq.(45):

∂(

ρ

h

)

∂

T

=

0

⇔

h0

/

cp

=

2T

−

ρ

1

/(

ρ

1

−

ρ

0

)

. (58)

The testsshow that the numerical scheme is stableif and only if thiscriterion is satisfied everywhere in the domain, regardlessoftheorderofthetime-steppingscheme,ortheorderofextrapolationforthepredictor.Furthermore,theresults showthattheminimalextrapolationordersinTable3needtobereachedinordertoachievethelowesterrors,buthigher ordersofextrapolationhavenoeffect.

7. Testcasesforthecoupledtransportequations

Thenumericalmethodforthecoupledtransportequationsandourcomputerimplementationareverifiedandvalidated withthefollowingtestcases.Wefirstinvestigateaseriesofmanufacturedsolutionsforthesystemoftransportequations, rangingfromconstantproperties(Section7.2), tovariabletransportproperties(Section7.3.1),andavariabledensity (Sec-tion7.3.2).AsinSection6,weleaveouttheunitsofthevariablesforthesemanufacturedsolutions.Section 7.4thenshows numericalresultsforbothconstant-densityandvariable-densityflowsthatcanbecomparedtoexperiments.

7.1. Implementation

Thesimulationsareperformedwithanin-housesolver

DGFlows

.Weuseahierarchicalsetofmodalbasisfunctions,so thatthereare



_P+P_d



degreesoffreedominad-dimensionalelementwithapolynomialorder

P

.

Allintegralsareevaluatedwithaquadraturesetthatissufficientlyaccuratetonegatethepolynomialaliasingeffectthat hasplaguedother DGsolvers.(See,e.g.,[27].)TheabscissaandtheweightsaretakenfromSolinet al. [44].Westore the valuesandderivativesofthebasisfunctionsonthequadraturesetforafastevaluationofintegralsandnumericalsolutions. Weverifiedthatnoneoftheresultsinthispaperchangedwhentheaccuracyofthequadraturewasincreased.

Fig. 2. Error(T−Tex_/Tex_{)att=1 forthetestcaseinSection6asafunctionoftheenthalpyoffseth}

0,usingmethod#1.Theverticaldottedlinesbound thevaluesforwhichEq.(57) heldatalltimesteps.

Allmeshesweregeneratedwiththeopen-sourcesoftwaretool

Gmsh

[16].ThemeshispartitionedwithMETIS[22]. WeusetheMPI-basedsoftwarelibrary

PETSc

tosolvethelinearequations[3,2].Thepressureequationissolvedwith a conjugate gradient method, whichis parallelized with an additive Schwarz preconditioner. The submatrix within each processispreconditionedwithanincompleteCholeskydecomposition.Thelinearsystemsfortheenthalpyandmomentum equationsaresolved witha GMRESmethod,whichisparallelized withan additiveSchwarzpreconditioner,and precondi-tionedwithanincompleteLUdecompositionwithineachprocess.

7.2. Taylor-GreenVortex

ThefirstmanufacturedsolutionforthetransportequationsistheTaylor-GreenVortex,whichisincompressibleandhas constant fluid properties.Wegeneralize thiswell-known analytical solutiontoinclude apassive scalartemperaturefield. Theenthalpyish

=

cpT .Theexactsolutionis

uex

=

exp



−

2˜t

 −

cos



˜

x



sin



y

˜



+

sin



˜

x



cos



y

˜





, pex

= −

ρ

4exp



−

4˜t

 

cos



2˜x



+

cos



2

˜

y



, Tex

=

exp



−

2˜t

/

Pr



cos



x

˜



cos



y

˜



,

(59)

onthedomainx

,

y

∈ [−

L

,

L

]

withDirichletboundaryconditionsand0

<

t

≤

1,where

˜

t

:=

ν

t

(

L

/(

n

π

))

2 ,

˜

x

:=

x L

/(

n

π

)

, y

˜

:=

y L

/(

n

π

)

, (60)

Fig. 3. Error(T−Tex_/Tex_{)att=1 forthetestcaseinSection6asafunctionoftheenthalpyoffseth}

0,usingmethod#2.Theverticaldottedlinesbound thevaluesforwhichthestabilitycriterioninEq.(58) isviolatedatsometimet.

Forthemanufacturedsolutionsinthissectionandthenext,thereportederrorsarerelativeintheL2-norm,thatis,



φ

− φ

ex



₂

φ

ex

_

2

=











T∈T



T

(φ

− φ

ex

)

2





(φ

ex

)

2 (61)

foraquantity

φ

.Eachintegralinthenumeratorisevaluatedwithanumericalquadrature,resultinginalargesumoverthe squaresofsmallnumbers.Wefoundthatanaiveimplementationgivesverylargeroundingerrorsthatdistorttheobtained ordersofconvergence.Wethereforeperformthedoublesummationovertheelementsandthequadraturepointswiththe Kahansummationalgorithm[48] anda128-bitfloatingpointnumber.

Fig. 4 shows the temporal convergence for L

=

1, n

=

1,

μ

=

0

.

01,

ρ

=

1, Pr

=

100, and a fourth-order polynomial spaceforthemassflux(i.e.,

P

m

=

4).WeperformedthesamenumericalexperimentsbyindependentlyvaryingthePrandtl

number(toPr

=

1), andthepolynomialorder(to

P

m

=

2),allofwhichyieldedsimilarresults.Allerrorssaturateatsmall

timesteps,wherethespatialdiscretizationerrordominates.

The Taylor-Greenvortex isofcoursea strangetest case,inthat itdoesnot featureanyscaleseparation:its frequency spectrum is comprised of delta Dirac functions. The solution is an eigenfunction of the diffusive terms, so there is no interactionbetweenthetransport terms.Furthermore,ithasDirichletboundaryconditionswherethereisoutflow,sothat thecontinuityequationisover-constrained.Thismanifestsitselfinastifflinearsystemforthepressure,thoughitsubiquity intheliteraturesuggeststhattheTaylor-Greenvortexisveryeasytosimulate.Thenextsectionfeaturesmorechallenging manufacturedsolutions.

7.3. Variable-propertymanufacturedsolutions

We presenttwo moremanufacturedsolutions withvariablefluid propertiesforthesystemofequations(1a)-(1c):one withaconstantdensity(inSection7.3.1),andonewithavariabledensity(inSection7.3.2).TherelationbetweenT andh

Fig. 4. Convergencetowardthe2DTaylorvortex(Eq.(59))attime˜t=1 withtemporalrefinementformesheswithN2_{squareelements.Theblackdashed} linesindicateidealsecond-orderconvergenceinδt.

The manufactured solutionisconstructed by workingbackward fromtheexact massflux andpressure,whichcan be chosen arbitrarily.Thechoiceofthepressureisoflittleconsequence,thoughwemake surethat bothm and p vary non-linearly in time, andthat they do not lie in the numerical solution space. The density then follows by integrating the continuityequationovertime.Thetemperatureandtheenthalpyareafunctionofthedensity.Theexternalforceandheat source cannowbe calculatedfromEqs.(1b) and (1c), whichwe dosymbolicallywiththePython SymPylibrary.We use polynomialsforthemanufacturedsolutionstokeepthesesymboliccalculationsfeasible.

Thisapproach differsslightlyfrompreviousliterature (e.g.,[43]),inthat wemake noattempttofindcleveranalytical solutions tothesystem(1) with avariabledensity,insteadsimplyintegratingthecontinuityequation toobtaina density thatconformstothemassflux.Somepreviouswork(e.g.,[39])hasincludedasourceterminthecontinuityequation,but thisappearslesssuitable foratime-splittingmethod,wherethecontinuityplays acentralroleinthediscretization(asin Section5.2),andwedonotwanttoadaptthenumericalschemetoconformwiththemanufacturedsolution.Ourapproach permitsarbitraryfunctionsfor

μ

=

μ

(

T

)

andk

=

k

(

T

)

.Weare notawareofprevious workwithamanufacturedsolution that is (i) compressible, (ii) satisfies the unmodified continuity equation, and(iii) use temperature-dependent transport properties.

Thedomainis

(

0

,

L

)

×(−

1

,

1

)

;seeFig.5.WeletL

=

10 inallcalculations,andusesquareelements.Theinflowboundary atx

=

0 hasDirichletboundaryconditions.Thegoalofthemanufactured solutionsisobviouslynot tomodelaparticular physicalphenomenon,butourconfigurationisreminiscentofapipe flowwithwallsat y

= ±

1 thatisheated asymmetri-cally,resultinginskeweddensityandvelocityprofiles.

All integrals are evaluated with a quadrature set with the usual polynomial accuracy of

(

3

P

m

−

1

)

, and we verified

that thisissufficienttointegrate uptomachine precisionby comparingtheresultswithahigher-orderquadraturesetof polynomialaccuracy

(

3

P

m

+

10

)

.

Fig. 5. Domain of the manufactured solutions in Section7.3.

Fig. 6. Constant-density manufactured solution in Eq. (62) at t=1.

We useaNeumannboundarycondition forthetemperatureattheoutlet,because thisisthe mostcommonchoice in practicalapplications.Theimposedheatflux(qN)followsfromtheknownexactsolution.WealsotriedimposingaDirichlet boundaryconditionforthetemperature,andfoundthatitmakesanegligibledifferenceinthenumericalerrors.

7.3.1. Constant-density,variable-propertymanufacturedsolution

Wefirstinvestigatethetemporalconvergenceforacasewherethetransportproperties(

μ

,k)dependonthe tempera-ture,butthedensityisconstant.Wechoosethepolynomialmanufacturedsolution

mex

= (

1

+

t3

)

1 L3

y



2x3

/

3

−

Lx2



(

y

−

1

)(

y

+

1

)

x

(

L

−

x

)



+

2 0



, pex

= (

1

+

t3

) (

L

−

x

)

3, Tex

= (1

+

t3

) (

2

−

y

)((

x

−

L

)/

L

)

2

/

6 (62)

with0

<

t

≤

1,whichsatisfies m2

=

0 on y

= ±

1, and

∂

kmk

=

0.The additionof theconstant

[

2

,

0

]

tomex ensures that

m1

>

0 everywhere,sothatthereisnobackflowattheoutlet, andnooutflowatanyoftheDirichletboundaryconditions. The densityis

ρ

=

1,but thetransport properties arenon-trivial:

μ

=

0

.

1

+

T

(

1

−

T

)

andk

=

μc

p

/

Pr,withPr

=

1.The

Fig. 7. Convergenceofthenumericalsolutiontowardtheconstant-densitymanufacturedsolution(Eq.(62))withtemporalrefinement.Thecharacteristic elementlengthisinverselyproportionaltoNy.Theblackdashedlinesindicateidealsecond-orderconvergenceinδt.

Table 4

Convergencetowardtheconstant-densitymanufacturedsolutioninEq.(62) (Fig.6)with spa-tialrefinementforthemixed-ordercase(Pm=2 andPh= Pp=1)withδt=2−12.

Ny Error T Conv. T Error u Conv. u Error p Conv. p

21 _{1.08e-2 - 1.11e-1 - 3.62e-3} -22 _{2.65e-3 2.03 1.83e-2 2.61 8.75e-4 2.05} 23 _{6.63e-4 2.00 2.62e-3 2.80 2.15e-4 2.02} 24 _{1.66e-4 2.00 3.52e-4 2.90 5.35e-5 2.01} 25 _{4.15e-5 2.00 4.57e-5 2.94 1.34e-5 2.00} 26 _{1.04e-5 2.00 6.43e-6 2.83 3.34e-6 2.00}

Fig.7displaysthetemporalconvergence.Weconsidervariousmeshes,varyingthenumberofelementsinthey-direction

(Ny), andthe spatial polynomial orders. The equal-order casedoes not appear to suffer fromthe inf-sup instability for

small

δ

t. Thevelocity andthetemperatureconvergewithsecond-orderaccuracyin

δ

t untiltheerrorsaturateswhen the spatialerrorstartstodominatethetemporalerror.Theorderofconvergenceforthepressureisslightlylower,intherange [1.5,2.0]. These orders of convergence for the velocity and pressure agree with what is found in previous literature on constant-propertyincompressibleflows(e.g.,[41],[15]).

ThespatialratesofconvergenceareinTable4.Asthemeshisrefined,allquantitieswithpolynomialorder

P

converge with

O





P+1



,where



∝

1

/

Ny isthe characteristicmeshlength. The convergencerateofthevelocity saturates athigh