HLLC-type methods for compressible two-phase flow in ducts with discontinuous area changes

Alexandra Metallinou Log

^a,∗

, Svend Tollak Munkejord

^b

, Morten Hammer

^b

aNorwegian University of Science and Technology, Department of Physics, Department of Energy and Process Engineering, Trondheim NO-7491, Norway

bSINTEF Energy Research, P.O. Box 4761 Torgarden, Trondheim NO-7465, Norway

a rt i c l e i nf o

Article history:

Received 18 December 2020 Revised 16 April 2021 Accepted 26 May 2021 Available online 29 May 2021 Keywords:

Finite-volume method HLLC solver Compressible flow Non-conservative system Two-phase flow Variable cross-section Nozzle flow

a b s t r a c t

Inthiswork,theHarten-Lax-vanLeerContact(HLLC)approximateRiemannsolverisextendedtotwo- phaseflowthroughductswithdiscontinuouscross-sections. Twomain strategiesareexploredregard- ingthetreatmentofthenon-conservativetermarisinginthegoverningequations.Inthefirst,labelled HLLC+S, the non-conservativeterm is discretized separately. In the second, labelled HLLCS, the non- conservativetermisincorporatedintheRiemannsolver.Themethodsareassessedbynumericaltests forsingleandtwo-phaseflowofCO₂,thelatteremployingahomogeneousequilibriummodelwherethe thermodynamicpropertiesarecalculatedusingthePeng–Robinsonequationofstate.Themethodshave differentstrengths,butingeneral,HLLCSisfoundtoworkbest.Inparticular,itisdemonstratedtobe equallyaccurateandmorerobustthanexistingmethodsfornon-resonantflow.Itisalsowell-balanced forsubsonicflowinthesensethatitconservessteady-stateflow.

© 2021 The Authors. Published by Elsevier Ltd.

ThisisanopenaccessarticleundertheCCBYlicense(http://creativecommons.org/licenses/by/4.0/)

1. Introduction

The simulation of two-phase flow through ductswith discon- tinuouscross-sectionsisessentialinseveralindustrialapplications.

Such simulationsareneededformodellinge.g. two-phaseflow in wellbores inthe oiland gasindustry [1],nuclear reactor coolant flows[2],emergencyventingofhydrocarbonpipelines[3]andcav- itation inrefrigeration systems[4].Systemslike thosementioned above can often be modelledasquasi one-dimensional withdis- continuous changes in cross-sectional area of the flow. The sys- temofequationsmodellingsuchflowcontainsanon-conservative term, andthis termcomplicatesnumericalsimulations greatly as itcancausenumericaloscillations[5,6]anddivergence[5].

Several authors have constructed numerical methods for the compressiblenozzleflowequations[1,5,7–9],andsystemsofsim- ilar form [10–16], developing “well-balanced” [17,18] schemes to capturetheflowbehaviouratdiscontinuities.Mostoftheearlyre- search has focused on thespecial caseof single-phaseflow with the ideal gas equation of state (EOS). Notable schemes include KrönerandThanh’swell-balancednumericalschemebasedonthe Lax-Friedrichsflux[19],whichwasextendedforresonantcasesin [5],Rochetteetal.’sVFRoebasedscheme[6]andCuongetal.’sGo-

∗Corresponding author.

E-mail addresses: [email protected] , [email protected] (A.M.

Log).

dunovschemebasedonan exactRiemannsolver[8].Brownetal.

[20] proposed thefirst methodologyfor resolvingtwo-phase CO₂ flow inpipeswithdiscontinuouscross-sectional areachanges for thehomogeneous equilibriumtwo-phase flow model(HEM)with the Peng-Robinson (PR) EOS [21] using the AUSM⁺-up scheme.

Recently, Abbasi etal. [1]developed a Godunov-type schemefor thetwo-phasedrift-fluxmodelwithvariablecross-section,though withsimpleEOSsforliquidandgas.

A HLLC-type method hasyet to be tested on the problem of compressible flow with discontinuous cross-sections. Note, how- ever,thattheHLLC-schemehasbeenextendedfortheEulerequa- tions in ducts ofsmoothly varying crosssections [22].HLLC-type schemesapplyinformationabouttheeigenstructureofthegovern- ingequationsintheirsolution[10,11,23],makingtheschemesless dissipativethangeneralmethodssuch asAUSM⁺-up [24].Forthe applicationontwo-phaseflow,theHLLC-scheme’saccurateresolu- tionofcontactdiscontinuities[23]isparticularlydesirableasthis alsomakesthe schememoreaccurateinresolving transitionsbe- tweengas,liquid,andmixtureflows. Astheeigenstructureofthe one-dimensional compressibleduct flow equationsis known, the advantagesabovemotivates theconstruction ofaHLLC solverfor thissystem.

Thisisfurther motivatedasaugmentedversionsof HLLChave beenconstructedforsimilarsystems,whereabruptchangesareac- counted for[10,11]. An augmented version of HLLC forthe Baer- Nunziato(BN)equations[25]wasdeveloped byTokareva andToro

https://doi.org/10.1016/j.compfluid.2021.105023

0045-7930/© 2021 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license ( http://creativecommons.org/licenses/by/4.0/ )

[10], giving promising results for many test cases. The method involves a nonlinear systemwhich was further linearizedby Lo- chonetal.[26].MurilloandGarcía-Navarro[11]alsodevelopedan augmented version ofHLLCforthe shallow-waterequations.This method produced promising results as well, though the authors note difficulties such as the need fora “source-fix” to avoid un- physicalsolutionsincertaincases.

Thecontributionofthisworkistodevelopandinvestigatetwo modifiedHLLCsolversforcompressibleductflowandassesstheir strengths and weaknesses. In particular, we will consider two- phase flow of CO₂, due to its use as a natural working fluid in refrigerationengineering, andtheimportance ofsafeandefficient CO2 transportation as part of CO2 capture and storage (CCS) as a climate-change mitigation technology [27]. We show that the presentmethodisbothrobustandaccuratewhensolvingchalleng- ingtwo-phaseRiemannproblems.

We willfirst presentthe equation system in more detail and briefly discuss theRiemann problemforthesystemin Section 2. The HEMandthe PR EOSareoutlined in Section3. Thenumeri- cal methodsarederivedinSection4,themethodsareassessedin Section5,andfinallysomeconcludingremarksandsuggestionsfor furtherworkaregiveninSection6.

2. GoverningequationsandtheRiemannproblem

The system of equations describing compressible one- dimensional flow of a single fluid in a rigid duct of variable cross-sectionalarea,A,is

Ut+F

(

^U

)

x=S, (1)

where

U=

⎛

⎜ ⎝ ρ

^A

ρ

^uA

EA A

⎞

⎟ ⎠

^{, F}

⁽

^U

⁾

⁼

⎛

⎜ ⎝ ρ

^uA

( ρ

^u2+p

)

^A

(

^E+p

)

^uA 0

⎞

⎟ ⎠

^{, S=}

⎛

⎜ ⎝

0 p^∂_∂^A_x

0 0

⎞

⎟ ⎠

^.

Here,

ρ

^{isthedensity,uthevelocity, E}=

ρ

(^e+¹₂u²)^thetotalen- ergy,ethespecificinternalenergy,andpthepressureofthefluid.

S is a non-conservative term. The set of Eq. (1) belongs to the classofnon-conservativeresonantsystems[5,28]meaningthatthe waveswhichariseinthissystemcaninteractand“resonate” with each other.Forsmooth solutions,thesystem(1)canberewritten inquasi-linearform,

Ut+A

(

^U

)

^Ux=0, (2)

whereAisthe Jacobianmatrixofthesystem. Notethatthenon- conservative termhas now been moved to the left-hand side of theequation.AfullderivationofAforageneralEOScanbefound in[29,AppendixD],andwehaveincludedthefullexpressionofA inAppendixA.

Itcanbeshown[19,28,29]thattheeigenvaluesofAare;

λ

^0=0,

λ

^1=u−c,

λ

^2=u,

λ

^3=u+c.

Notethatanyoftheeigenvalues

λ

1,

λ

2,

λ

3maycoincidewith

λ

0, giving rise toresonance inthe system[28].The system ofequa- tionsishyperbolicawayfromthepointswhere

λ

1=

λ

0or

λ

3=

λ

0

andnonstrictlyhyperbolicwhen

λ

2=

λ

0[28]. 2.1. TheRiemannproblem

ConsidertheRiemannproblemforcompressibleductflow,

Ut+F

(

^U

)

^x=S, (3)

U

(

x,0

)

=

UL, ifx<0

UR, ifx≥0 , (4)

where U_L and U_R are two different constant states. A thorough analysis on the characteristic fields, Riemann invariants and the solutionto thisRiemann problemispresented by Andrianovand Warneckein[28].

WhenthereisnochangeinA,A_L=A_R,thesystem(1)reduces totheEulerequations.Wethenhavethesamecharacteristicsand RiemanninvariantsasfortheEulerequationsassociatedwiththe eigenvalues

λ

1,

λ

2,

λ

3.TheRiemanninvariantsare

s, u+2c

^across

dx1

dt =u−c (5)

u,p across dx₂

dt =u (6)

s,u−2c

^across

dx3

dt =u+c, (7)

wheresisthespecificentropyand^{isthefirstGrüneisenparam-} eter,

⁼

ρ

¹

∂

^p

∂

^e

ρ= 1

ρ

^cv

∂

^p

∂

^T

ρ. (8)

Here,c_visthespecificheatcapacityatconstantvolume.Admissi- blewaves forthesolution tothe Riemannproblemare then rar- efactionsandshocksassociatedwith

λ

1,

λ

3 anda contactdiscon- tinuityassociatedwith

λ

2.

Atpointswithdiscontinuousarea change,thereisastationary contactdiscontinuityassociatedwiththeeigenvalue

λ

0=0[19],the 0-wave.Acrossthe0-wavewehavethefollowingRiemanninvari- antsasshownin[28]

A

ρ

^{u, s, h+1}₂^{u2, across d}_d^x_t^{0 =}

λ

^{0=0, (9)}

where h=e+_ρ^p is thespecific enthalpy of thefluid. The invari- ants describe the conservationof massflux, entropy andstagna- tion enthalpy overthe area change. Theaddition of thiswave in thesolutiontotheRiemannproblemcausescomplicationssuchas non-uniqueness[28] andresonance[5,19,30].InFig.1weprovide anexampleofthestructureofaRiemannproblemsolutioninthe caseof subsonic flow i.e.

|

^u

|

<c fromleft to right. The example wascreatedusingAndrianov’sprogram[31](CONSTRUCT).

3. Thermodynamicmodels

Inthiswork,wemodelthefluidasanidealgasforbenchmark testsof numericalsolvers of theequation system(1),defined by theequationofstate(EOS)

p=

ρ ( γ

⁻¹

)

e, (10)

where

γ

^{istheratioofspecificheats}

γ

= ^cc^pv.Inadditiontobench- marktesting, it is alsorelevant to studythe system(1)for two- phase flowof liquidandgas.To modelthis, we apply thehomo- geneousequilibriummodel(HEM)withthePeng-Robinson(PR)EOS [21].ThePREOSisgivenby

p= RT

v

m−b−

α

^a

v

²m+2b

v

m−b², (11) where

v

misthespecificmolarvolumeofthefluidandRisthegas constant.a,band

α

^{aredefinedas}

a=0.45724R²T_c²

pc , (12)

b=0.07780RTc

pc, (13)

Fig. 1. The characteristics of a Riemann problem for 1D compressible duct flow with subsonic flow where ρL> ρR, p L> p Rand A L> A R, giving a rarefaction to the left, a stationary contact discontinuity (red), then a contact discontinuity and a shock to the right. Created using [31] . (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 2. A one dimensional interval subdivided into grid cells, j, with cell centers at x jand faces x _j−1/2, x _j+1/2.

and

α

=

1+

0.37464+1.54226

ω

−0.26992

ω

²

1−

T Tc

2

,

(14) where Tc, pc and

ω

^{arethe critical} temperature, criticalpressure andtheacentricfactorofthespecies.ForCO₂,theseare

pc=7.3773MPa, Tc=304.35K, and

ω

=0.2236. (15) ThePREOSonlygivesresidualheatcapacities,c^res_p ,c^res_v .Inorderto compute the total heat capacities cp=c^ideal_p +c^res_p ,c_v in JK⁻¹kg⁻¹ weusethefollowingestimate,

c^ideal_p =479.107+1.524318T−1.078176·10⁻³T²+

+3.38976·10⁻⁷ T³+2.8876·10⁻¹¹ T⁴. (16) In theHEMitis assumedthatthe twophasesare inthermal, chemical andmechanicalequilibrium,whichisvalidifthephases are well-mixed. Mixture properties are then used in the flow Eq. (1). Inthiswork, SINTEF’s thermodynamic library[32,33] has been applied to providesolutions for the HEMwiththe PR EOS.

Detailsonthespecificmethodsappliedinthelibrarytoobtainrel- evant variables are presentedin [34], though we ignoreherethe presenceofanysolid.

4. Numericalmethods

The computational domain isdiscretized in finitevolumes j

asdepictedinFig.2.Weusetwo differentkinds offinite-volume methods(FVMs)tosolveEq.(1)onthisgrid.ThefirstFVMisanal- ogous tothespatialdiscretizationthatBrownetal.applyin[20], withanEulertimestepgiving

Uⁿ_j⁺¹=Uⁿ_j −

^t

^x

Fj+1/2−Fj−1/2

+

^tSj, (17)

Fig. 3. The flux function F ⁺_j−1/2approximates the flux F _j−1⁺ _/2just to the right of the interface at x _j−1/2. The flux function F ⁻_j+1/2approximates the flux F _j+1⁻_/2just to the left of the interface at x j+1/2.

where F=F(^UL,U_R) ^{is a numerical flux function} approximating theaverage fluxF atthecell interfacesx=x_j−1/2,x=x_j+1/2, and S_j approximatesthe contributionofthenon-conservative termin cell j.

The second FVMis aconservative Godunovscheme whichin- cludes thenon-conservative termin thenumericalflux functions [11].TheFVMtakesthefollowingform

U_jⁿ⁺¹=Uⁿ_j −

^t

^x

F⁻_j+1/2−F⁺_j−1/2

, (18)

whereagainanEulertimestepisusedforthetemporaldiscretiza- tion. Here, F^±=F^±(^UL,U_R,S) ^{are numerical flux functions ap-} proximatingtheaverageflux, F, rightnextto theeast, F_j⁻_+1/2,and west,F_j−1⁺ _/2,cellfacesasillustratedinFig.3.

Inthe following,we will briefly review theHLLC methodand thensuggesttwomodifiedHLLC-typemethodstoapproximatethe fluxesforthecompressibleductflow.

4.1. TheHLLCapproximateRiemannsolver

TheHLLC method,proposedby Toro,SpruceandSpeares [23], approximatesthecellinterfaceRiemannproblembyathree-wave solution;

U

(

^x/t

)

=

⎧ ⎪

⎨

⎪ ⎩

UL, ifx<

v

Lt, U_L^HLLC, if

v

Lt≤x<

v

Ct, U_R^HLLC, if

v

Ct≤x<

v

Rt, UR, ifx≥

v

Rt,

(19)

where

v

L and

v

R are thefastestsignal velocities arising fromthe initial conditionofthe Riemann problem, and

v

C is thespeed of thecontactwave.TheintermediatestatesU_L^HLLC,U_R^HLLC areapprox- imatedtobeconstant,

U_L^HLLC= 1

^t

( v

C−

v

L

)

tvC

^tvL

U

(

x,

^t

)

^{dx (20)}

and

U_R^HLLC= 1

^t

( v

R−

v

C

)

tvR

tvC

U

(

x,

^t

)

^dx, (21)

theyarehoweverunknownandmustbeestimated.HLLCapproxi- matesthenumericalfluxfunctionby

Fj+1/2=

⎧ ⎪

⎨

⎪ ⎩

FL, if0<

v

L, F_L^HLLC, if

v

L≤0<

v

C, F_R^HLLC, if

v

C≤0<

v

R, FR, if0≥

v

R.

(22)

Theintermediatestatefluxes,F_L^HLLCforpositivesubsonicflow,and F_R^HLLC for negative subsonic flow, are also unknown. In order to determine the fluxes, Rankine-Hugoniot (RH) relations are used across the waves and the additional set of Riemann invariants acrossthecontactdiscontinuityisappliedtoclosethesystem.The RHrelationstatesthatacrossawave

^F=

v

^U, (23)

where

v

isthespeedofthewave.Forcompressibleduct flow,we findthroughsomemanipulationthattheintermediatefluxesF_K^HLLC, K=L,Rcanbeexpressedas

F_K^HLLC=FK+

v

K

(

^UK^HLLC−UK

)

, (24) wheretheintermediatestatesareapproximatedby

U_K^HLLC=

ρ

KAK

v

K−uK

v

K−

v

C

¹

v

C EK

ρK+

( v

C−uK

)

v

C+_ρK(v^p_K^K_−uK)

, (25)

K=R,L, (25)

and

v

C= p_R−p_L+

ρ

Lu_L

( v

L−u_L

)

−

ρ

Ru_R

( v

R−u_R

)

ρ

L

( v

L−u_L

)

⁻

ρ

R

( v

R−u_R

)

^{. (26)}

4.2. Wave-speedestimates

The HLLC solver needs estimates for the wave speeds

v

L and

v

R.Thereare severaldifferentapproachestoestimate thesewave speeds, some of which are outlined in [35], Section 10.5. In this work, the Roe average wave speed estimate [36] is used. Both Davis[37]andEinfeldt[38]suggestusingtheRoeaveragedeigen- valuesforthewavespeeds;

v

L,j+1/2=min

λ

1

(

^Uj

)

,

λ

1

(

^Uj+1/2

)

,

v

R,j+1/2=max

λ

³

(

^Uj+1

)

,

λ

³

(

^Uj+1/2

)

, (27)

whereUistheRoeaverageoftheconservedvariables.TheRoeav- eragedvariablescanbefoundbytheRoeaveragedmatrixA(^UL,U_R) [36],whichmustsatisfycertainconditions.

We follow the approachof Evje andFlåtten [39] andMunke- jord [40] for the two-fluid model, which also involves a non- conservative term, andsearch fora Roeaveraged matrixAwhich satisfiesthefollowingconditions:

R1A(^UL,U_R)(^UR−U_L)=^F(^UL,U_R)

R2A(^UL,UR)^hasrealeigenvaluesandisdiagonalizable,and R3A(^UL,U_R)→A(^U)^smoothlyasUL,U_R→U,

wherein^F(^UL,U_R)^{isformulatedas}

^F

(

^UL,UR

)

=

⎛

⎜ ⎝

{ ρ

^uA

} { ( ρ

^u2+p

)

^A

}

^−pˆ

{

^A

}

{ (

^E+p

)

^uA

}

0

⎞

⎟ ⎠

^{. (28)}

Here,

{

^x

}

^=xR−x_L, (29)

and pˆ is a particular average of the pressures from the left and rightstatespˆ=pˆ(^UL,U_R)^,similarlyto

α

k(^UL,U_R)^in[39],[40].

Acan be determined by finding a specialaverage of thestate vectorsU_L andU_R,U(^UL,U_R)^{,such thatA}=A(^U)^{, pˆ}=pˆ(^U)^.Aset ofaveragessatisfyingR1–R3are:

ρ

^A=

ρ

^LAL+

ρ

^RAR

2 , (30)

Aˆ=AL+AR

2 , (31)

ˆ u=

ρ

LA_Lu_L+

ρ

RA_Ru_R

ρ

LA_L+

ρ

RA_R , (32)

Hˆ=

ρ

^LALHL+

ρ

^RARHR

ρ

^LAL+

ρ

^{RAR , (33)}

whereH_k=h_k+¹2u²_k,k=L,R.

4.3. HLLCwithaddednon-conservativeterm,HLLC+S

TheHLLC schemeassumesathreewave solution,howeverwe can still apply the scheme to compressible duct flow provided that we also account for the fourth, stationary wave. We apply the FVM (17) with the HLLC numerical flux function. This FVM requires a representation of the non-conservative term, S_j. The discretizationof thistermrequiresspecialcare toensure numer- ical stability. We follow the approach of Brown et al. [20] for their AUSM⁺-up scheme and apply a discretization of the non- conservativetermwhichsatisfiesthenon-disturbancerelationdis- cussed by Liou et al. [41]. The relation states that under steady conditionswithu=0andp=const.

∂ (

^Ap

)

∂

^{x =p}

∂

^A

∂

^{x. (34)}

The following discretization, which satisfies the non-disturbance relation,isused:

Sj= pj

^x

⎛

⎜⎝ 0 Aj−A_j−1

0 0

⎞

⎟⎠^{, ifuj>0andSj=pjx}

⎛

⎜⎝ 0 A_j+1−Aj

0 0

⎞

⎟⎠^,ifuj≤0.

(35) 4.4. HLLCSapproximateRiemannsolver

We willhere derive an augmented version ofHLLC, following in part the approach of Murillo and García-Navarro [11] for the shallow-water equations and the approach of Tokareva and Toro [10]forthe Baer-Nunziatoequations.Wefollowthe namingcon- vention of Murillo and García-Navarro [11] and call this method

“HLLCS”,emphasizingthatourmethodisverysimilartotheirdis- cretizationofthesourcetermintheshallow-waterequations.The HLLCS approximateRiemannsolverassumesa four-wave solution instead of a three-wave solution, incorporating the 0-wave. Sim- ilarly to HLLC, we will assume that the waves separate constant intermediatestates.

ForsystemsintheformofEq.(1),MurilloandGarcía-Navarro derivedthefollowingconsistencyconditionwhichtheapproximate intermediatestatesmustsatisfy:

1

^t

( v

R−

v

L

)

_tvR tvL

U

(

^x,

^t

)

^dx=

v

RUR−

v

LUL−

(

^FR−FL

)

+S

v

R−

v

L , (36)

Fig. 4. Integration control volume [ x L, x R] ×[0 , t] in the x −t plane. The control volume contains the two fastest signal velocities, vLand vRfrom the Riemann prob- lem. The solution consists of three inner states separated by the stationary wave at x = 0 and the contact discontinuity of positive speed, vC.

where S= 1

^t xR

xL

t 0

Sdt dx. (37)

Two differentestimates of S are used in thiswork and they are presentedinSection4.4.4.HLLCSwillbedevelopedtoensurethe subsoniccasesatisfiesthecondition(36).Forsupersonicflow,the fluxesareeasilyfoundaswillbeshownbelow.

4.4.1. Supersonicflow

Forpositivesupersonicflow,theflux justtotheleft ofthein- terface,x=0,issimplyF_L,giving

F⁻_j+1/2=F_L. (38) The flow just to the right of the interface has passed the area changesuchthat

F⁺_j+1/2=FL+S. (39) Similarly fornegativesupersonicflow, thenumericalfluxesat^t become:

F⁻_j+1/2=FR−S, (40)

F⁺_j+1/2=FR. (41) 4.4.2. Subsonicflow

An illustrationofa controlvolume containingthewave struc- tureofaRiemannproblemforpositivesubsonicflowisshownin Fig. 4. In this casethere are three unknown intermediate states separatedbythestationarywaveatx=0andthecontactdiscon- tinuity,U_L⁻,U_R⁺andU_R⁺⁺.

Weapproximatetheintermediatestates,U_L⁻,U_R⁺andU_R⁺⁺by U_L⁻=₋¹_tvL

₀ tvL

U

(

x,

^t

)

^dx U_R⁺=¹_tvC

tvC 0

U

(

x,

^t

)

^dx U_R⁺⁺=_t(v¹_R−vC)

tvR

^tvC

U

(

x,

^t

)

^dx

⎫ ⎪

⎪ ⎪

⎪ ⎪

⎬

⎪ ⎪

⎪ ⎪

⎪ ⎭

. (42)

In ordertoestimate theintermediate fluxes, theRH condition is applied across all thewaves in the problem. The RH relations are

F_L⁻−F_L=

v

L

(

^UL⁻−U_L

)

^{, (43)}

Fig. 5. Integration control volume [ x L, x R] ×[0 , t] in the x −t plane. The control volume contains the two fastest signal velocities vL, vRfrom the Riemann problem.

The solution consists of three inner states separated by the stationary wave at x = 0 and the contact discontinuity of negative speed, vC.

F_R⁺−F_L⁻−S=

v (

^UR⁺−U_L⁻

)

=0, (44)

F_R⁺⁺−F_R⁺=

v

C

(

^UR⁺⁺−U_R⁺

)

^{, (45)}

FR−F_R⁺⁺=

v

R

(

^UR−U_R⁺⁺

)

. (46) ItcanbeshownthattheRHrelations(43)-(46)areenoughtosat- isfytheconsistencycondition(36).Toclosethesystem,weimpose theRiemanninvariantsacrossthestationarywaveandthecontact discontinuity,

u⁺⁺_R =u⁺_R =

v

C

p⁺⁺_R =p⁺_R

, (47)

(

^A

ρ

^u

)

⁻L =

(

^A

ρ

^u

)

⁺R

s⁻_L =s⁺_R

(

^u₂²+h

)

⁻L =

(

^u₂²+h

)

⁺R

}

. (48) TheRHconditionacrossthewaveassociatedwiththewave speed

v

L gives

ρ

L⁻=

ρ

L

v

L−uL

v

L−u⁻_L , (49)

p⁻_L =pL+

ρ

^L

( v

L−uL

)(

^u−_L −uL

)

, (50)

E⁻_L =

ρ

^L

v

L−uL

v

L−u⁻_L

_E

L

ρ

^L+

⁽

^{u−L −uL}

⁾

^{u−L +}

pL

ρ

^L

( v

L−uL

)

, (51)

andthe RH conditionacross the wave associated withthe wave speed

v

R gives

ρ

R⁺⁺=

ρ

^R

v

R−u_R

v

R−u⁺⁺_R , (52)

p⁺⁺_R =pR+

ρ

R

( v

R−uR

)(

^u++R −uR

)

, (53)

E⁺⁺_R =

ρ

R

v

R−uR

v

R−u⁺⁺_R

_E

ρ

RR+

(

^u++R −uR

)

^u++R + pR

ρ

R

( v

R−uR

)

.

(54)

Fig. 6. The graphs show for which values of p ⁻_L, p ⁺_R that f 1= 0 (blue) and f 2= 0 (red) for Mod. A (a) and Mod. B (b). For Mod B, an estimate of the point where |^f| is minimized is marked with an x. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Eqs.(49)-(51)and(52)-(54)withtheRiemanninvariantsconstitute a nonlinear systemwhich canbe solved iteratively.Tokareva and Toro[10]obtainedasimilar, butlargersystemofequationswhich mustbesolvedfortheBaer-Nunziatoequations.

Both forcompressible duct flow andthe Baer-Nunziato equa- tions,eitherthepressures p⁻_L, p⁺_R =p⁺⁺_R orthevelocitiesu⁻_L,u⁺_R = u⁺⁺_R can bechosen as independent variables to solvethe system.

As statedin [10],there is no difference between the approaches fromatheoreticalpointofview asthetworepresentationsofthe system are mathematically equivalent. Following Tokareva et al.

[10], we choose p⁻_L,p⁺_R as the independent variables to ensure pressurepositivitywhensearchingforsolutionsofthesystem.

Wethereforeexpressu⁻_L andu⁺_R =u⁺⁺_R using p⁻_L andp⁺_R, u⁻_L

(

^p−L

)

=uL+ p⁻_L −pL

ρ

L

( v

L−uL

)

^{, (55)} u⁺_R

(

^p+R

)

=uR+ p⁺_R−pR

ρ

R

( v

R−uR

)

^{. (56)}

We then havethatU_L⁻=U_L⁻(^p−L)^{, suchthat s−}L =s⁻_L(^p−L)^,anden- forcing the Riemann invariant s⁻_L =s⁺_R =s, we have that s⁺_R = s⁺_R(^p−L)^{.The relationformassflux andthe relationforstagnation} enthalpythengivethefollowing:

f=

A_Lρ_L⁻(^p−_L)^u−_L(^p−_L)−A_Rρ_R⁺ p⁺_R,s(^p−_L)

u⁺_R(^p+_R) h⁺_R

p⁺_R,s(^p−L) +¹2

u⁺_R(^p+R)2

−! h⁻_L

p⁻_L,s(^p−L) +¹2

u⁻_L(^p−L)2"

=0.

(57) Thesearetwoequationsforthetwoindependentvariablesp⁻_L,p⁺_R. Thesystem(57)canbesolvediterativelybye.g.Newton-Raphson’s methodanditmayhavezeroorupto threesolutions.Ifthe sys- temhasmultiplesolutions,wechoosethesolutionwhichsatisfies thefollowingcriteria:

C1 The solution is self-consistent in the sense that the Riemann problemforthestatesU_L⁻(^p−L),U_R⁺(^p−L,p⁺_R)^{providewavespeed} estimateswhichsuggestsubsonicflow.

C2 The solution has the highest entropy s⁻_L(^p−_L)=s⁺_R =s of the self-consistentsolutions.

If there areno solutions, we approximate p⁻_L,p⁺_R as thepoint whichminimizestheabsolutevalue f₁(^p−L,p⁺_R)+f₂(^p−L,p⁺_R)^{of f.}

Oncep⁻_L andp⁺_R aredetermined,thestateU_L⁻canbecalculated using Eqs. (55), (49) and (51). Withthis we can finally find the unknownfluxesF_L⁻andF_R⁺ fromEq.(43)andEq.(44),giving F_L⁻=FL+

v

L

(

^UL⁻−UL

)

^{, (58)}

F_R⁺=F_L⁻+S. (59)

Fig. 7. Comparison of the exact density solution (black line) and the solutions of HLLC+S (red circles), HLLCS with RS (blue plus signs) and HLLCS with FS (green crosses) at t/t ref= 0 . 02 for Test 1. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Thenegativesubsonicflowcasecanbeseenassimplythemir- rorimage ofthepositiveflowcase.We nowhavethestatesU_L⁻⁻, U_L⁻andU_R⁺as illustrated inFig. 5.

Anequivalentsystemto(57)canbefoundforthiscaseandthe samecriteriaC1andC2canbeappliedtochooseavalidsolution.

4.4.3. Solutionforstationarywaves

Supposenow thatwe havethestatesU_L,U_R whichsatisfy the conditionsforastationarywaveacrosstheareachange,

(

^A

ρ

^u

)

L=

(

^A

ρ

^u

)

R, hL+u²_L

2 =hR+u²_R

2, sL=sR. (60) TheexactsolutionfortheRiemannproblem(3)–(4)withthetwo statesU_L,U_RisajumpfromU_L toU_Rattheareachange.Thesolu- tionwhichsatisfiesthecriteriaC1andC2isp⁻_L =p_L andp⁺_R =p_R. TheintermediatestatesthenbecomeU_L⁻=UL andU_R⁺=U_R⁺⁺=UR. This means that forstationary waves, when the correct solution ischosen,theintermediatestatesfoundintheHLLCSapproximate Riemannsolverareexact.

4.4.4. Thenon-conservativetermforHLLCS

Inthiswork,twonon-conservativetermsaretestedtoestimate the fluxesbased on the HLLCS approximateRiemann solver. The

Fig. 8. Results of the convergence test for HLLC+S (red line with circles) and HLLCS with RS (blue line with plus signs) for Test 1. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 9. Density solution of HLLC+S (red dashed line), HLLCS with RS (blue dotted line) and HLLCS with FS (green dash-dotted line) compared to the exact solution (black line) for Test 2 at t/t ref= 0 . 1 , with N cells= 10 0 0 . (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

firstapproximatenon-conservativetermisgivenby

SRS=

⎛

⎜ ⎝

0 ˆ

p

(

^AR−AL

)

0 0

⎞

⎟ ⎠

^{, (61)}

where pˆ is the Roe-averaged pressure introduced in Section 4.2. We thereforecallthistheRoe-average-basedterm(RS). RSis for- mulated generallysuch that itmaybe appliedon subsonic,sonic andsupersonicflow.

Forsubsonicflow,thenonlinearsystemofequationsdetermin- ing theapproximate intermediatestatesandfluxesis solved.The non-conservative term is then given implicitly by the RH condi- tions.Forpositivesubsonicflow, wegetthatthenon-conservative termmustbe

S_FS+=F_R⁺⁺−

v

C

(

^UR⁺⁺−U_R⁺

)

^−FL⁻. (62) Similarlyfornegativesubsonicflow,wegetthat

SFS-=F_R⁺−F_L⁻⁻+

v

C

(

^U_L⁻−U_L⁻⁻

)

. (63) Asthenon-conservativetermincludestheapproximatefluxes,we callittheflux-basedterm(FS).Notethatthisestimate onlyholds iftheHLLCSapproximateRiemannsolverhasasolution.FSisonly formulated for subsonic flow andmay therefore only be applied forsubsonicflowproblems.

4.4.5. TheHLLCS-basedfluxes

TheHLLCSmethodapproximatesthefluxfunctionsF⁻_j+1/2and F⁺_j+1/2 neededfortheFVM(18)asshowninAlgorithm1.

Algorithm1:TheHLLCS solver.Ifsubsonic flowisidentified, asolveriscalledtofindavalidsolutionsatisfyingC1andC2 or an optimization method is used to minimize f. When a solutionisfound,

v

C andtheintermediate statesU_L⁻,U_R⁺and U_L⁻⁻orU_R⁺⁺are returned.

Result:FluxesfortheHLLCSsolver,F_R⁺andF_L⁻. if

v

L>0then

F_L⁻=F_L F_R⁺=F_L⁻+S end

if

v

L≤0and

v

R>0then

callsolver,returning

v

Candintermediatestates;

if

v

C≥0then

F_L⁻=F_L+

v

L(^U_L⁻−U_L) F_R⁺=F_L⁻+S

else

F_R⁺=F_R−

v

R(^UR−U_R⁺) F_L⁻=F_R⁺−S

end end

if

v

R≤0then F_R⁺=F_R F_L⁻=F_R⁺−S end

SetF⁻_j+1/2=F_L⁻andF⁺_j+1/2=F_R⁺.

Remark: Note that for a (subsonic) steady-state wave across thearea change,applyingS_FS willgive that F⁻_j+1/2=F_L⁻=FL and F⁺_j+1/2=F_R⁺=F_R.InsertingthisintheFVM(18),wefindthat

U_jⁿ⁺¹=Uⁿ_j

∀

j, (64)

i.e.theHLLCS-basedFVMwithFSconservesthesteady-statesolu- tionexactly.ThismeansthattheFVMiswell-balanced.

4.5. Summary

Inthiswork,weapplytwo finite-volumemethodsHLLC+Sand HLLCS.ThenumericalschemefortheHLLC+Ssolverisgivenby U_jⁿ⁺¹=Uⁿ_j −

^t

^x

Fj+1/2−Fj−1/2

+

^tSj, (65)