A novel pressure-free two-fluid model for one-dimensional incompressible multiphase flow

Delft University of Technology

A novel pressure-free two-fluid model for one-dimensional incompressible multiphase flow

Sanderse, B.; Buist, J. F.H.; Henkes, R. A.W.M.

DOI

10.1016/j.jcp.2020.109919

Publication date

2021

Document Version

Final published version

Published in

Journal of Computational Physics

Citation (APA)

Sanderse, B., Buist, J. F. H., & Henkes, R. A. W. M. (2021). A novel pressure-free two-fluid model for

one-dimensional incompressible multiphase flow. Journal of Computational Physics, 426, 1-18. [109919].

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

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.

Contents lists available atScienceDirect

Journal

of

Computational

Physics

www.elsevier.com/locate/jcp

A

novel

pressure-free

two-fluid

model

for

one-dimensional

incompressible

multiphase

flow

B. Sanderse

a

,

∗

,

J.F.H. Buist

a

,

c

,

R.A.W.M. Henkes

b

,

c a_{CentrumWiskunde&Informatica(CWI),Amsterdam,theNetherlands}

b_{ShellTechnologyCentreAmsterdam,Amsterdam,theNetherlands} c_{DelftUniversityofTechnology,Delft,theNetherlands}

a

r

t

i

c

l

e

i

n

f

o

a

b

s

t

r

a

c

t

Articlehistory: Received14April2020

Receivedinrevisedform30August2020 Accepted9October2020

Availableonline22October2020 Keywords:

Two-fluidmodel Pressure-freemodel Constraint Runge-Kuttamethod

A novel pressure-free two-fluid model formulation is proposed for the simulation of one-dimensional incompressible multiphase flow in pipelines and channels. The model is obtainedby simultaneouslyeliminatingthe volume constraintand the pressure from thewidelyusedtwo-fluidmodel(TFM).Theresulting‘pressure-freetwo-fluidmodel’ (PF-TFM)hasanumberofattractivefeatures:(i)itfeaturesfourevolutionequations(without additional constraints) that can be solved very quickly with explicit time integration methods; (ii) it keeps the conservation properties ofthe original two-fluid model,and thereforethecorrectshockrelationsincaseofdiscontinuities;(iii)itssolutionssatisfythe twoTFMconstraintsexactly:thevolumeconstraintandthevolumetricflowconstraint;(iv) itoffersaconvenientformtoanalyticallyanalysetheequationsystem,sincetheconstraint hasbeenremoved.

Astaggered-gridspatialdiscretizationandanexplicitRunge-Kuttatimeintegrationmethod areproposed,whichkeeptheconstraintsexactlysatisfiedwhennumericallysolvingthe PF-TFM.Furthermore,forthecaseofstronglyimposedboundaryconditions,anoveladapted Runge-Kuttaformulationisproposedthatkeeps thevolumetricflowexactintimewhile retaining high order accuracy. Severaltest cases confirm the theoretical properties and showtheefficiencyofthenewpressure-freemodel.

©2020TheAuthor(s).PublishedbyElsevierInc.Thisisanopenaccessarticleunderthe CCBYlicense(http://creativecommons.org/licenses/by/4.0/).

1. Introduction

The two-fluidmodel(TFM)forone-dimensionalmultiphase flowisan importantmodeltostudy,forexample,the be-haviour of oiland gastransport in long pipelines. Inthis articlewe study its incompressiblesimplification. Anexample application is flushing (cleaning) an oil-filled pipeline at atmospheric conditions with water, in which case the incom-pressible assumption isaccurate. Like insingle-phaseflow models,the incompressibility assumptionnecessitates the use ofcarefulspace- andtime-discretization methods.IntheincompressibleTFM,the maindifficultyindefining accurate nu-mericsarisesfromthepresenceofthevolumeconstraint(thetwo phasesshould togetherfillthepipeline),its associated volumetricflowconstraint(themixturevelocityfieldshouldbedivergence-free),andtheroleofthepressure.

OneapproachfornumericallysolvingtheincompressibleTFMistoeliminatethepressurefromthefour-equationsystem (+1constraint)andtorewritethissystemintoatwo-equationsystem(+2constraints).Thisleadstothe‘no-pressurewave’

*

Correspondingauthor.

E-mailaddress:B.Sanderse@cwi.nl(B. Sanderse). https://doi.org/10.1016/j.jcp.2020.109919

0021-9991/©2020TheAuthor(s).PublishedbyElsevierInc.ThisisanopenaccessarticleundertheCCBYlicense (http://creativecommons.org/licenses/by/4.0/).

modelorthe‘fixed-flux’model(FFM)suggestedbyWatson[22],andusedforexamplein[1,6,12,14].Anothertwo-equation modelisthe‘reverseddensity’modeldevelopedbyKeyfitzetal.[10] andemployedforexamplein[21].Inthesemodels, thepressure isgenerallycomputedasapost-processing step.Ageneralproblemwiththesetwo-equationsystemsisthat theyareonlyvalidforsmoothsolutions.Thismeansthatinthepresenceofshocksdifferentjumpconditionsareobtained thanforthefour-equationmodel[1] (anexampleofshockwavesinthetwo-fluidmodelarerollwaves[22]).Furthermore, thefixed-fluxassumptionoftenlimitsthesestudiestostationaryboundaryconditions.

Anotherapproachis tostick tothefour-equation formulationandkeepthepressure, andto useapressure-correction method.Apressureequationisthentypicallyobtainedbysubstitutingthemomentumequationsinthecombinedmass con-servationequation,eitherwithemployingthevolumeconstraintequation[11],orwithoutemployingthevolumeconstraint equation[8,12,20].Inourrecentwork[17] weproposedaconstraint-consistentpressureequation,whichwasshowntobe bothageneralizedRiemanninvariantofthecontinuousequations,andahiddenconstraintofthesemi-discretedifferential algebraicequation system.TheproposedRunge-Kuttatime integrationmethodwasshowntobehigh-orderaccuratewhile satisfyingboth thevolumeconstraintandthevolumetricflowconstraint. However,thisapproachiscomputationally rela-tivelycostly,asoneneedstosolveapressurePoissonequationateachstageoftheRunge-Kuttamethod,withtheadditional difficultythatthePoissonmatrixdependsontheactualsolutionandthereforechangesintime.

In thispaperwe propose a novel incompressibletwo-fluid model formulationinwhich the pressureand thevolume constrainthavebeeneliminated.The eliminationprocess willbesuch thatthe resultingpressure-freemodelstillsatisfies boththevolumeconstraintandthevolumetricflowconstraint.Comparedtothefour-equationapproach,ournovel formu-lationismuchfastertoevaluatesinceitdoesnotrequirethesolutionofaPoissonequation.Comparedtothetwo-equation approach, the newmodel keepsthe correctshock solutions, andfurthermore itdoes not rely onthe assumption of the ‘fixed-flux’.

The outline of this paperis as follows:in section 2 the incompressibletwo-fluid model equations are given, in sec-tion3thenewpressure-freemodelisderived,andinsections4and5appropriate spatialandtemporaldiscretizationsare proposed.Resultsareshowninsection6.

2. Governingequationsforincompressibletwo-phaseflowinpipelines

The incompressible two-fluid modelcan be derived by considering the stratifiedflow of liquid andgas ina pipeline (for arecent discussionofthe two-fluid model,seeforexample [12]).The mainassumptions that wemake are thatthe flowisone-dimensional,stratified,incompressible,andisothermal.Thetransversepressurevariationisintroducedvialevel gradient terms.Surface tensionis neglected.This leads tothe followingsystemofequations,which we callthe ‘original’ TFM[17]:

∂

U

∂

t

+

∂

f

(

U

)

∂

s

+

h

(

U

)

∂

p

∂

s

+

S

(

U

)

=

0

,

(1) where U

=

⎛

⎜

⎜

⎝

U1 U2 U3 U4

⎞

⎟

⎟

⎠ =

⎛

⎜

⎜

⎝

ρg

Ag

ρ

lAl

ρg

ugAg

ρl

ulAl

⎞

⎟

⎟

⎠ ,

f

(

U

)

=

⎛

⎜

⎜

⎝

ρ

gugAg

ρ

lulAl

ρg

u2 gAg

+

Kg

ρl

u2 lAl

+

Kl

⎞

⎟

⎟

⎠ ,

h

(

U

)

=

⎛

⎜

⎜

⎝

0 0 Ag Al

⎞

⎟

⎟

⎠ ,

S

(

U

)

=

⎛

⎜

⎜

⎝

0 0 Sg Sl

⎞

⎟

⎟

⎠ ,

(2)

supplementedwiththevolumeconstraintequation:

Ag

+

Al

−

A

=

0

.

(3)

In these equations the subscript denotes either gas (g) or liquid (l). The model features four evolution equations, one constraintequation,andfiveunknowns(U and p),whichareafunction oftheindependentvariabless (coordinatealong thepipelineaxis)andt (time).Theprimitivevariables Ag,ug etc.canbeexpressedfullyintermsofU ,e.g.

Ag

=

U1

ρ

g

,

ug

=

U3 U1

.

(4)

ρ

denotesthedensity(assumedconstant), A thecross-sectionalareaofthepipe, Ag and Al (alsoreferred toasthe hold-ups) thecross-sectionalareasoccupiedby thegasorliquid, R the piperadius,h the height oftheliquidlayermeasured fromthebottomofthepipe,u thephasevelocity, p thepressureattheinterface,

τ

theshearstress (withthewallorat the interface), g the gravitationalconstant,

ϕ

thelocalupward inclination ofthepipeline withrespectto thehorizontal, andthecomponentsofgravityare gn

=

g cos

ϕ

andgs

=

g sin

ϕ

.SeeFig.1.

ThelevelgradienttermsK forapipe(compressibleorincompressibleflow)aregivenby[16]:

Kg

= −

ρ

ggn



(

R

−

h

)

Ag

+

1 12P 3 gl

,

(5) Kl

= −

ρ

lgn



(

R

−

h

)

Al

−

1 12P 3 gl

,

(6)

Fig. 1. Stratified flow in a pipeline. Left: cross-sectional view; right: side view.

whileforachannel(compressibleorincompressible)onehas: Kg

= −

1 2

ρ

ggnA 2 g

,

(7) Kl

=

1 2

ρ

lgnA 2 l

.

(8)

In apipe geometry, Al andh arerelatedby a non-linearalgebraic expression,see e.g. [17];ina channelone has Al

=

h (assumingunitdepth).Thewettedandinterfacialperimeters Pg,PlandPgl canbeexpressedintermsoftheliquidhold-up

Al ortheinterfaceheighth (see[17]).

Thesourceterms Sg andSl constitutefriction,gravity,anddrivingpressuregradient,givenby Sg

=

τ

glPgl

+

τg

Pg

+

ρg

Aggs

+

Ag dp0 ds

,

(9) Sl

= −

τ

glPgl

+

τ

lPl

+

ρ

lAlgs

+

Al dp0 ds

.

(10)

Thefrictionmodelsaredescribedin[17].ThesourcetermsSgand Sl donotcontainspatialortemporalderivativesofthe unknownsofthesystem(U and p).Thedrivingpressuregradient dp0

ds isimposedincaseofperiodicboundaryconditions. Initialandboundaryconditionsdeterminewhethertheequationsformawell-posedinitialboundaryvalueproblem.Itis well-knownthatthisisnotalwaysthecase:thetwo-fluidmodelequationscanbecomeill-poseddependingonthevelocity differencebetweenthephases[3,13,18].Inthispaperwewillrestrictourselvestotestcasesforwhichthemodelis well-posed.Furthermore,aswill becomeclearin thenext section, werequirethe initial conditionstobe consistent, meaning thatatt

=

0 thesolutionshouldnotonlysatisfythevolumeconstraintbutalsothevolumetricflowconstraint:

∂

∂

s

ugAg

+

ulAl



=

0

.

(11)

3. Anewpressure-freetwo-fluidmodel

3.1. ThehiddenconstraintsoftheincompressibleTFM

In thissection we show how thepressure andthe constraintcan be removedfromthe incompressibleTFM.For this purpose,we use thecharacteristic analysisofthe incompressibletwo-fluid model proposed in[17]. In that work,it was shownthatthevolumeconstraint,equation(3),couldberewrittenintwoalternativeforms.

Thefirstformisthevolumetricflowconstraint,obtainedbypremultiplyingtheTFMequationswith

l3

=



1 ρg 1 ρl 0 0

,

(12) leadingto l3

·



∂

U

∂

t

+

∂

f

(

U

)

∂

s

+

h

(

U

)

∂

p

∂

s

+

S

(

U

)



=

0

,

⇒

∂

Ag

∂

t

+

∂

Al

∂

t

+

∂

∂

s

Agug

+

Alul



=

0

,

⇒

∂

∂

s

Agug

+

Alul



=

0

.

(13)

The last equationis thevolumetricflowconstraint. Notethat inthisderivationwe used thetime-derivativeofthevolume constraint,

∂

∂

t

Ag

+

Al



=

0

,

(14)

sothevolumetricflowconstraintcanbeinterpretedasacombinationofmassequationsandvolumeconstraint.Integrating thevolumetricflowconstraintinspacegives

Agug

+

Alul

=

Q

(

t

).

(15)

Thesecondformisthepressureequation,obtainedbypremultiplyingtheTFMequationswith

l4

=



0 0 _ρ1 g 1 ρl

,

(16) leadingto l4

·



∂

U

∂

t

+

∂

f

(

U

)

∂

s

+

h

(

U

)

∂

p

∂

s

+

S

(

U

)



=

0

,

⇒

∂

∂

t

ugAg

+

ulAl



+

∂

∂

s



u2_gAg

+

Kg

ρ

g

+

u_l2Al

+

Kl

ρ

l



+



Al

ρ

l

+

Ag

ρg



∂

p

∂

s

+

Sg

ρg

+

Sl

ρ

l

=

0

,

⇒ −



Al

ρl

+

Ag

ρg



∂

p

∂

s

= ˙

Q

(

t

)

+

∂

∂

s



u2_gAg

+

Kg

ρg

+

u 2 lAl

+

Kl

ρl



+

Sg

ρg

+

Sl

ρl

.

(17) Attheheartofthisderivationisthetime-derivativeofequation(15),whichinitselfwasobtainedfromthetime-derivative ofthevolumeconstraint,seeequation(14).Inotherwords,equation(17) canbeinterpretedasanequationforthe second-ordertime-derivativeoftheconstraint:

∂

2

∂

t2

Ag

+

Al



=

0

.

(18)

Notethatthepressuregradientequation (17) isdifferentfromtheexpressioncommonlyusedinliterature(seee.g.[7]),in whichthegasandliquidmomentumequationsareaddedwithoutscalingbytheirrespectivedensities.Althoughsuchapressure equation can be usefulfor postprocessingcomputations,it cannot beused toreplace thevolume constraint, becausethe volumeconstrainthasnotbeenemployedinitsderivation.Incontrast,ourderivedpressureequation(17) isacombination ofthemomentumequations,massequations,and thevolumeconstraint.

3.2. Thepressure-freeTFM

Thekeyinsightofthispaperisthefollowing.Theexpressionforthepressuregradient,equation(17),canbesubstituted backintotheTFMequations.Forthispurpose,firstwritethepressureequationinshortnotationas

∂

p

∂

s

= −

ρlρg

ˆ

ρ

l4

·



∂

f

(

U

)

∂

s

+

S

(

U

)



−

ρlρg

ˆ

ρ

Q

˙

(

t

),

(19) where

ˆ

ρ

(

U

)

:=

ρ

gAl

+

ρ

lAg

.

(20)

ThepressuregradienttermaspresentintheTFMthenreads

h

(

U

)

∂

p

∂

s

=

⎛

⎜

⎜

⎝

0 0 Ag Al

⎞

⎟

⎟

⎠

∂

_∂

p_s

= −

⎛

⎜

⎜

⎜

⎝

0 0 Agρlρg ˆ ρ Alρlρg ˆ ρ

⎞

⎟

⎟

⎟

⎠

·

l4

·



∂

f

(

U

)

∂

s

+

S

(

U

)



−

⎛

⎜

⎜

⎜

⎝

0 0 Agρlρg ˆ ρ Alρlρg ˆ ρ

⎞

⎟

⎟

⎟

⎠

Q

˙

(

t

)

= −

B

(

U

)

·



∂

f

(

U

)

∂

s

+

S

(

U

)



−

⎛

⎜

⎜

⎜

⎝

0 0 Agρlρg ˆ ρ Alρlρg ˆ ρ

⎞

⎟

⎟

⎟

⎠

Q

˙

(

t

),

(21) with B

(

U

)

=

⎛

⎜

⎜

⎜

⎝

0 0 Agρlρg ˆ ρ Alρlρg ˆ ρ

⎞

⎟

⎟

⎟

⎠

·

l4

=

1

ˆ

ρ

⎛

⎜

⎜

⎝

0 0 0 0 0 0 0 0 0 0 Agρl Agρg 0 0 Al

ρ

l Al

ρ

g

⎞

⎟

⎟

⎠ .

(22)

TheimportantconclusionisthatthepressuregradienttermintheTFMcanbewrittenintermsofalinearcombinationoftheflux

vectorandthesourceterms.Uponsubstitutingtheexpressionforh

(

U

)

∂_∂p_s intotheTFM,weobtainthenewpressure-freeTFM

(PF-TFM)as PF-TFM:

∂

U

∂

t

+

A

(

U

)

∂

f

(

U

)

∂

s

+

A

(

U

)

S

(

U

)

+

c

(

U

) ˙

Q

(

t

)

=

0

,

(23) where A

(

U

)

=

I

−

B

(

U

)

=

⎛

⎜

⎜

⎝

1 0 0 0 0 1 0 0 0 0 1

−

Agρl

/

ρ

ˆ

−

Ag

ρg

/

ρ

ˆ

0 0

−

Alρl

/

ρ

ˆ

1

−

Alρg

/

ρ

ˆ

⎞

⎟

⎟

⎠ ,

c

(

U

)

=

⎛

⎜

⎜

⎜

⎝

0 0

−

Agρlρg ˆ ρ

−

Alρlρg ˆ ρ

⎞

⎟

⎟

⎟

⎠

.

(24)

Equation (23) isa closedsystemoffourequations infourunknowns(providedthat Q

˙

(

t

)

isknown):boththe constraint and the pressure havebeen removed. The volume constraint is ‘builtinto’ thissystem ofequations (it was used inthe derivationofthepressureequation),anddoesnotneedtobeprovidedasafifthequation.Wenotethatmakingthemodel pressure-free comesata small‘price’:oneneeds Q

˙

(

t

)

asan input tothesystemofequations.In other words:removing the pressurerequires theprescription of Q

˙

(

t

)

.Fortunately, Q

˙

(

t

)

equals zeroinmanycases, e.g. incaseofsteadyinflow conditionswithprescribedliquidandgasflowrates.Forunsteadyinflow,e.g.anincreasingliquidandgasproduction, Q

˙

(

t

)

is knownfromtheboundary condition specification.Only inspecialcases,such asperiodicboundary conditions, Q

˙

(

t

)

is unknown;we will assume itto be zeroin thatcase andinvestigateby comparingPF-TFM solutions toTFM solutionsto whatextentthisisindeedtrue.

3.3. PropertiesofthePF-TFM

ThePF-TFMsharesthepropertyoftheTFMthatitcannotbewritteninfullconservativeform.IntheTFMthismanifests itself in thepresence ofthe hold-up fractionsin front ofthe pressure gradientterm; in thePF-TFM thisenters through thepresenceofthe A

(

U

)

matrix.Notethat A

(

U

)

isasingularmatrix,whichimpliesthatthesystemcannotbesimplified by pre-multiplyingwiththeinverseof A.The eigenvaluesof A

(

U

)

∂_∂f(_UU) are exactlythesameasthoseobtainedfromthe analysis ofthe original TFM inprimitive variables in[17], see Appendix A; asa consequence,the PF-TFMhas the same characteristicvelocitiesastheTFM.

ThePF-TFMdoesnotrequirethesolutionofaPoissonequationforthepressure,whichsignificantlyacceleratesthecode (this will be shownin the results section). Givena solution U tothe PF-TFM, the pressure can be obtainedas a post-processingstepbysolvingequation(19).AnexistingTFMimplementation(basedonapressurePoissonequation)caneasily beadaptedbysimplyscalingthecomputedfluxes∂_∂f_s usingtheweightsfromthe A matrix,andaddingtheQ

˙

(

t

)

term.The pressurePoissonequationcanthenbeskipped.

Comparedtothefixed-fluxmodel(FFM),thereare twoimportantadvantages. First,incontrasttotheFFM,thePF-TFM retainsthe correctshockrelations (i.e.,thoseofthe TFM).Second, the FFMis typicallyused with‘fixed-flux’assumption

Q

(

t

)

=

constant [12],whereasthisassumptionisnotmadeinourmodel.

Although the constraint is ‘built into’ the pressure-free TFM, upon discretization care should be taken, because the volume constraint is not built into the PF-TFM in its original form, equation (3), but in a differentiated sense, namely equation(18).Thiswillbediscussednext.

4. Spatialdiscretizationonastaggeredgridtopreserveconstraints

Inthederivationofthepressure-freemodel,equation(23),thefollowingimportantpropertywassilentlyused:theterm presentinthetime-derivativeinthemomentumequations,e.g.U3

=

ρ

gAgug,isthesametermasinthefluxterminthe massequations.Thispropertyneeds tobesatisfiedatadiscrete levelinordertoobtainapressure-freeTFMthat satisfies thevolumeconstraint;thisisachievedherebytheuseofastaggeredgrid,seeFig.2.

Thestaggered gridconsistsofNp

=

N ‘pressure’ andNu

=

N

+

1 ‘velocity’volumes(withslightchangesdependingon theboundaryconditions).Althoughthepressureisnotpresentinthenewformulation,thisconceptissowell-knownthat westick totheterminology.Themidpointsofthevelocityvolumeslie onthefacesofthepressurevolumes.Thepressure, hold-up andphase masses are definedin the centreof the pressure volumes, whereas the velocity andmomentum are definedin the centreofthe velocity volumes. Fordetails we refer to [16]. Theunknowns are thevector of conservative variables U

(

t

)

∈ R

2Np+2Nu_{,whichincludethefinitevolumesizes(weusethenotation U todistinguishfromthevectorU}

thatappearsinthecontinuousequations):

U

(

t

)

=

⎛

⎜

⎜

⎝

U₁

(

t

)

U₂

(

t

)

U₃

(

t

)

U₄

(

t

)

⎞

⎟

⎟

⎠ =

⎛

⎜

⎜

⎝

[(

ρg

Ag



s

)

1

. . . (

ρg

Ag



s

)

N

]

T

[(

ρ

lAl



s

)

1

. . . (

ρ

lAl



s

)

N

]

T

[(

ρ

gAgug



s

)

1/2

. . . (

ρg

Agug



s

)

N+1/2

]

T

[(

ρl

Alul



s

)

1/2

. . . (

ρl

Alul



s

)

N+1/2

]

T

⎞

⎟

⎟

⎠ .

(25)

Fig. 2. Staggered grid layout including left boundary.

Includingthevolumesizeinthevectorofunknownsleadstoaclearphysicalinterpretation(e.g.U₁ hasunitsofmass)and allowsgeneralizationtothecaseoftime-dependentfinite-volumesizes.U₁

,

U₂

∈ R

Np _andU

3

,

U4

∈ R

Nu containthevalues ofmassandmomentumatthemassandvelocityvolumes,respectively,andareonlyafunctionoftime.Withthesevectors, theexpressionforthegashold-upandgasvelocityis(similarto(4)):

A_g

(

t

)

=



−1 p U1

(

t

)

ρg

,

ug

(

t

)

=



−_u1U₃

(

t

)

Ip



−p1U1

(

t

)

,

(26)

where



p

∈ R

Np×Np isadiagonalmatrixthatcontainsthepressurevolumesizes,



u

∈ R

Nu×Nu adiagonalmatrix contain-ingthevelocityvolumesizes,andIp

∈ R

Nu×Np aninterpolationmatrixfrompressuretovelocitypoints.

We start withconservation of mass for the gas phase. Integration of the first equation in (23) in s-directionover a pressurevolumegives:

d dt

U1,i



+

U3,i+1/2



u,i+1/2

−

U3,i−1/2



u,i−1/2

=

0

.

(27)

Acrucialdetailinthisequationisthattheconvectivefluxesaredirectlyexpressedintermsofthemomentum U₃.Forall volumes,thediscretizationiswrittenas

dU₁

dt

=

F1

(

U

)

= −

Dp



−_u1U₃



,

(28)

where Dp isaspatialdifferencingmatrix,see[17].Forclarity:F1 denotesthefirstcomponentofthespatial discretization of A∂_∂f_s

+

A S

+

cQ .

˙

Forconservationofmomentumweproceedinasimilarway.Integrationofthethirdequationin(23) ins-directionover avelocityvolumegives:

d dt

U3,i+1/2



+

A33,i+1/2

(

f3,i+1

−

f3,i

)

+

A34,i+1/2

(

f4,i+1

−

f4,i

)

+

A33,i+1/2S3,i+1/2

+

A34,i+1/2S4,i+1/2

−



Agρl

ρg

ˆ

ρ



i+1/2



si+1/2Q

˙

(

t

)

=

0

,

(29) where the notation Ai j indicates an entry inthe A matrix. Note that this discretization wouldalso havebeen obtained whenfirstdiscretizingtheTFMandthenapplyingthepressure-freetransformationtoit.

Inthetestcasesinthisarticle,eitherupwind(first-orderaccurate)orcentral(second-orderaccurate)schemesareused toexpress f intermsoftheconservativevariables.Ingeneral,theschemereads

f3,i

= (

ρg

Ag

)

iu2g,i

+

Kg,i

.

(30)

In thecentral scheme we haveug,i

=

1₂

(

ug,i−1/2

+

ug,i+1/2

)

,whereas in the first-orderupwindscheme ug,i

=

ug,i−1/2 if

ug,i−1/2

>

0. Ag and ug follow fromequation (26), and Kg from equation (5). A similar expression holds for f4, with subscripts g replaced by l. For more details, we refer to [17]. The particular choice for f3 and f4 does not affect the derivationofthediscretepressure-freemodel.Ifdesired,higher-ordermethodsforapproximating f3 and f4 canbeeasily employedinstead.

Insummary,thediscretizationofthemomentumequationofthegasphaseiswrittenas dU₃

andsimilarforU₄.We stressthat incontrastto theoriginalTFM, F3 doesnotcontainthepressure,butinsteadcontains the A-matrixandvolumetricsourceterminvolving Q .

˙

Theentirespatial discretization(mass andmomentumequations) canbesummarizedas

dU

(

t

)

dt

=

F

(

U

(

t

),

t

).

(32)

Wewillnowcheckwhetherthevolumeandvolumetricflowconstraintsaresatisfiedwiththisspatialdiscretization.The first-ordertime-derivativeofthevolumeconstraintisgivenby

d dt

A_g

+

A_l



=

d dt





−_p1U₁

ρg

+



−_p1U₂

ρl



,

= 

−1 p d dt



U₁

ρg

+

U₂

ρl



,

= −

−1 p



F₁

ρg

+

F₂

ρ

l



,

= −

−1 p Dp





−_u1U₃

ρg

+



−_u1U₄

ρ

l



.

(33)

Thesecond-ordertime-derivativethenfollowsas d2 dt2

A_g

+

A_l



= −

d dt





−_p1Dp





−_u1U₃

ρg

+



−_u1U₄

ρ

l



,

= −

−1 p Dp



−u1 d dt



U₃

ρg

+

U₄

ρl



,

= −

−1 p Dp



−u1



F₃

ρg

+

F₄

ρ

l



.

(34)

Thelasttermcanbesimplifiedbyrealizingthatthe A matrixinthePF-TFMhasbeenconstructedsuchthat: F₃

ρg

+

F₄

ρl

= 

uQ

˙

(

t

),

(35)

where Q

˙

∈ R

Nu

=

_1Q

˙

₍

_t

₎

_{isavector with Q}

₍

_t

₎

_{ineachentry.Thecorrectnessofthisexpressioncanbe easilycheckedby}

substitutingtheexpressionsforF₃andF₄(seeequation(29)),andfor A andc,seeequation(24).Seealsoequation(B.12). Consequently,equation(34) gives

d2 dt2

A_g

+

A_l



= −

_p−1DpQ

˙

(

t

),

=

0

.

(36)

Thelastequalityholdsbecauseapplicationofthedifferencingmatrix Dp toaconstantvectoryields0.

Integratingequation (36) twiceintime,givenaninitialconditionthatsatisfiesthevolumeconstraint( A_g

(

t0

)

+

Al

(

t0

)

=

A),thenleadsto

A_g

(

t

)

+

Al

(

t

)

=

A

.

(37)

Inotherwords,thevolumeconstraintremainssatisfiedintimebyemployingaspatialdiscretizationonastaggeredgrid. 5. Explicittimediscretization:efficientandconstraint-consistent

5.1. Basicformulation

Whereas theprevious section showedthat the pressure-free TFM posedrequirements on the spatial discretization in orderto beconstraint-consistent, nofurtherrequirementsare necessaryforthe temporaldiscretization:a simpleforward Euler time discretization already keeps the constraint property. The extension to generic explicit Runge-Kutta methods is then straightforward. Note that due to the absence of the pressure, the spatially discretized system can be advanced straightforwardlyintime(withouttheneedtosolveaPoissonequation).

TheforwardEulerstepisgivenby

Animportanttermpresentin F

(

Un

)

isthediscreteapproximationtothelast terminequation (29),beingthevolumetric flowsourceterm.Wesplit F

(

U

)

as

F

(

U

)

= ˆ

F

(

U

)

+

c

(

U

) ˙

Q

(

t

).

(39)

TheforwardEulerdiscretizationofthelastterm,asgivenbyequation(38),yields:

option A:



tc

(

Un

) ˙

Q

(

tn

).

(40)

Analternativeformulation(optionB)ofthesourcetermwillbediscussedinsection5.2.

As before, the key question is whether the constraintis satisfied at the new time level tn+1 (assuming that it was satisfiedatthecurrenttimeleveltn_):

An_g+1

+

An_l+1

−

A

=



−1 p Un1+1

ρg

+



−_p1Un₂+1

ρl

−

A

,

= 

−1 p



Un₁

+ 

t F₁

(

Un

)

ρg

+

Un₂

+ 

t F₂

(

Un

)

ρl



−

A

,

=

An_g

+

An_l

−

A







0

−

t



−_p1Dp





−_u1Un 3

ρ

g

+



−_u1Un₄

ρ

l



,

= −

t



−_p1Dp



−u1



_Un 3

ρ

g

+

Un₄

ρl



.

(41)

Inordertosimplifythisequationfurther,weinserttheexpressionforthemomentumequationandemployproperty(35)

Un₃

ρg

+

Un₄

ρl

=

Un₃−1

+ 

t F₃

(

Un−1

)

ρ

g

+

Un₄−1

+ 

t F₄

(

Un−1

)

ρl

,

=

U n−1 3

ρg

+

Un₄−1

ρ

l

+ 

uQ

˙

(

tn−1

).

(42)

Here we usedoption Aforthediscretization ofthe volumetricsourceterm. Givenan initial condition forU₃ and U₄ at

t

=

t0 whichsatisfiesthevolumetricflowconstraint, Dp



−u1



U1₃

ρg

+

U1 4

ρ

l



=

0

,

(43)

andgiventhat DpQ

˙

(

t

)

=

0,thevolumetricflowconstraintcanbewrittenasatelescopicsumoverallprevioustimesteps:

Dp



−u1



Un₃

ρg

+

Un₄

ρl



= . . . =

Dp



−u1



U1₃

ρg

+

U1₄

ρl



=

0

.

(44)

Thevolumeconstraintatthenewtimelevel,equation(41),canthusbeevaluatedas

An_g+1

+

A_ln+1

−

A

=

0

.

(45)

Insummary,whenapplyingtheforwardEulermethodtoequation(32) andstartingfromconsistentinitialconditions,both thevolumeconstraintandthevolumetricflowconstraintaresatisfiedateachtimestep.

Note that neitherequation (42), norequation (44), guarantees that thesolution Un+1 will satisfy theexact flow rate

Q

(

tn+1

)

,sinceonlythetime-derivative Q

˙

(

tn

)

hasbeenused.However,when insteadofoptionA,we applythefollowing first-orderdiscretizationofthevolumetricflowsourceterm

option B: c

(

Un

)(

Q

(

tn+1

)

−

Q

(

tn

))

(46)

intoequation (42), not onlythevolumetricflow constraintissatisfied,butinadditionthe actualvalue ofthevolumetric flowstaysequaltothespecifiedvolumetricflowwhenmarchingintime:

Un 3

ρg

+

Un 4

ρ

l

= 

uQ

(

tn

).

(47)

Inotherwords,whenadaptingtheforwardEulerdiscretizationforthevolumetricflowsourceterm,itispossibletosatisfy the volumetric flow constraint and the actual value of Q

(

t

)

exactly. The generalization of this idea to high-order time integrationmethodsrequirescarefulconsiderations,aswillbediscussednext.

5.2. High-orderaccuracyandintegrationofQ

˙

(

t

)

Thehigh-ordertimeintegrationmethodsweconsiderareexplicitRunge-Kutta(RK)methods.ExplicitRKmethodsform anexcellent combinationofaccuracyandstabilityforconvection-dominatedsystemsand,incontrasttothestandardTFM, theycanbeappliedtotheproposedPF-TFMinastraightforwardmanner,becausetheconstrainthasbeeneliminated.When consideringoptionAforhandlingthevolumetricflowsourceterm,theresultingmethodreads:

option A Un,i

=

Un

+ 

t i−1



j=1 ai jF

(

Un,j

,

tj

),

Un+1

=

Un

+ 

t s



i=1 biF

(

Un,i

,

ti

),

(48) (49)

wherea andb arethecoefficientsoftheButcher tableau,and s thenumberofstages.Since explicitRKmethods canbe seenasacombinationofforwardEulersteps, itisstraightforwardtoprovethatthey satisfythetwo constraints,bothfor thestagevaluesUn,iandforUn+1.

When consideringoptionB forhandlingthevolumetricsourceterm, thereisan open questionatwhichtime levelto evaluatec

(

U

)

attheintermediatestages.Onepossibilitythatleadstoanexactvolumetricflowistotakeateachstage

c

(

Un

)(

Q

(

tn,i

)

−

Q

(

tn

)).

(50)

However, this choice doesnot fit inthe Runge-Kuttaframework, and theresulting method suffersfrom lossof order of accuracy.

The question how to employ optionB in conjunctionwitha high-order methodcan be answered with thefollowing insight.ThetimediscretizationcorrespondingtooptionB,equation(46),isobtainedwhenderivingthepressure-freemodelonthe

fullydiscretelevelinsteadofonthecontinuouslevel. In other words: option A isobtained by first eliminating the pressure

fromthetwo-fluidmodel,andthendiscretizingtheresultingsystem;optionBisobtainedbyfirstdiscretizingthetwo-fluid model,andtheneliminatingthepressure.ThedetailsofthisderivationareshowninAppendixB.Withthisinsight,a high-orderRKmethodforoptionBisstraightforwardtoconstruct:eliminatethepressurefromthe(fullydiscrete)high-orderRK methodsproposedin[17],andrewriteasatimeintegrationmethodforthePF-TFM.

Wefocusonthecases

=

3,i.e.athree-stagethird-ordermethod,forwhichthefollowingmethodresults:

option B Un,1

=

Un

,

Un,2

=

Un

+ 

t 1



j=1



a2 jF

ˆ

(

Un,j

,

tj

)



+

c

(

Un,1

)(

Qn,2

−

Qn

),

Un,3

=

Un

+ 

t 2



j=1



a3 jF

ˆ

(

Un,j

,

tj

)



+

c

(

Un,1

)

a31 a21

(

Qn,2

−

Qn

)

+

c

(

Un,2

)



(

Qn,3

−

Qn

)

−

a31 a21

(

Q n,2

₋

_Qn

₎



,

Un+1

=

Un

+ 

t 3



i=1



biF

ˆ

(

Un,i

,

ti

)



+

c

(

Un,1

)

b1 a21

(

Q n,2

₋

_Qn

₎

₊

c

(

Un,2

)



b2 a32

(

Qn,3

−

Qn

)

−

b2 a21 a31 a32

(

Qn,2

−

Qn

)



+

c

(

Un,3

)



(

Qn+1

−

Qn

)

−

b2 a32

(

Qn,3

−

Qn

)

−



b1 a21

−

b2 a21 a31 a32



(

Qn,2

−

Qn

)



.

(51) (52) (53) (54)

Itcanbecheckedthroughsubstitutionthatthevolumetricflowremainsequaltothespecifiedvolumetricflow Q

(

t

)

,both atthestagelevelsandatthenewtimelevel.Itisimportanttoremarkthat,althoughthisisstrictlyspeakingnota Runge-Kuttamethodappliedto(32),itis aRunge-KuttamethodappliedtotheoriginalTFM.Theorderofaccuracyofthismethod wasanalyzedin[17],andaspecificthird-ordermethodwasproposed:

c a b

=

0 0 1 2 1 2 1

−

1 2 1 6 2 3 1 6 (55)

This methodis such that the lower diagonal entries (a21, a32) are nonzero (thisis required inequations (53)-(54)),and order reductionthat can resultfromunsteadyboundary conditionsisavoided (detailsinnext section).The c-coefficients

are usedtodeterminetheintermediate timelevels (ti

₌

_tn

₊

_c

i



t)andshouldnotbe confusedwiththe vectorc or c as presentinthePF-TFM.

Notethatthetemporalaccuracy(thirdorder)ofthismethodishigherthanofthespatialdiscretization(secondorderfor thecentral scheme,andfirstorderfortheupwindscheme).Thisisbecausefirst- andsecond-orderRunge-Kuttamethods are notstableforpure convectionproblems(discretizedwithcentralschemes), independentofthe timestep,incontrast tothird-ordermethods.Furthermore,thethird-ordermethodleavesroom forfutureimprovementsintermsofhigh-order spatialdiscretizationmethods.

Althoughequations(51) - (54) might seemcumbersomefromanimplementationpointofview,one canusean exist-ing Runge-Kuttaimplementationandsimplyadd theadditionalvectorsonthe right-handside that involvec and Q .The effectivecoefficientsofthescheme,whichinvolvecombinationsofRunge-Kuttacoefficients,canbeeasilyprecomputed.

Insummary,wehaveproposedhigh-ordertimeintegrationofthePF-TFMforboth‘weak’(optionA)and‘strong’(option B)impositionofthevolumetricflowsourceterm.Inbothcasesthetwoconstraintsaresatisfied,andinthelattercasealso thevolumetricflowremainsequaltothespecifiedvalue.ThesecondapproachrequiresaspecializedRunge-Kuttamethod, which isslightlymore involvedto implement. Notethat in thespecial casethat Q

˙

(

t

)

=

0, we have F

ˆ

=

F andoptionB becomes a classic Runge-Kutta method,equivalent to option A. In thiscasewe willresort to a classic four-stage fourth orderRunge-Kuttamethod.

5.3. Boundaryconditions

The boundary condition treatment follows thecharacteristic approach outlined in[17],which has theadvantage that it doesnot requirethe pressurein theformulation. Here we shortlysummarizethe (important)caseof unsteadyinflow conditions,withprescribedgasandliquid(mass)flows,Ig

(

t

)

andIl

(

t

)

,respectively.Therearetwochoicesfortheboundary conditions:prescriptioninstrongform(prescribetheactualvalueofIgandIl ateachtimestep),orinweakform(prescribe thetime-derivatives

˙

Ig and

˙

Il andintegrateintimewiththeRunge-Kuttamethod).

It isimportanttonote that theboundaryconditionprescriptionshouldbeconsistentwiththediscretizationofthevolumetric

sourceterm c

(

U

) ˙

Q

(

t

)

in orderto make surethat Q remains uniformin space, so that the volumetricflow constraint is

satisfied.TheweakboundaryconditionprescriptionisconsistentwithoptionA.Thestrongboundaryconditionprescription isconsistentwithoptionB.

5.4. Keepingtheconstraintexactintimebypreventingmachineerroraccumulation

Inderivingequation(41) theerrorinthevolumeconstraintattheprevioustime stepwas imposedtobeexactlyequal to zero.In actual computations,thevalue of An

g

+

Anl

−

A might not exactly equal zerodueto thepresence ofmachine

precisionerrors.Bysimplyleavingthisterminthenumericalalgorithm,onecanavoidaccumulationofmachineprecision errorsthatcouldspoiltheaccuracyofthecomputations[5,19].Inotherwords,oneshouldnot imposethevolumeconstraint termin(41) tobezero,butinsteadsimplycomputeitsvaluebasedonthesolutionatthelasttimestep( An_g

+

An_l

−

A)and usethiswhiletimestepping.Inpracticeweachievethisby computingtheerrorinthevolumeconstraintaftereachstage oftheRunge-Kuttamethod:

ε

n_A,i

= 

−1 p



Un₁,i

ρg

+

Un₂,i

ρ

l



−

A

,

(56)

andaddingthistermtothesolutionaftereachstage:

ˆ

Un,i

=

Un,i

+

ε

n,i

,

(57) where

ε

n,i

_{= −}

1 2

⎛

⎜

⎜

⎝

ρg



p

ε

n_A,i

ρ

l



p

ε

n_A,i 0 0

⎞

⎟

⎟

⎠ .

(58)

Table 1

Parameter values for the test case withtheKelvin-Helmholtzinstability.

Parameter Value Unit

ρl 1000 kg/m3 ρg 1.1614 kg/m3 R 0.039 m g 9.8 m/s2 μg 1.8·10−5 Pa s μl 8.9·10−4 Pa s 10−8 _m L 1 m 6. Results 6.1. Kelvin-Helmholtzinstability

Weperformtheclassicviscous Kelvin-Helmholtzinstabilitytestcase,seee.g.[11,17] andTable1forparametervalues. We choosevaluesfor Al andul,ul

=

1 m

/

s and

α

l

=

Al

/

A

=

0

.

9,andthen computethe gasvelocity andthebody force necessarytosustainthesteadysolution:

ug

=

8

.

01 m

/

s

,

dp0

ds

= −

87

.

87 Pa

/

m

.

(59)

Thevelocitydifference ug

−

ul isbelowthe‘Kelvin-Helmholtz’instabilitylimit[11,17],whichmeansthat,atleastinitially, theinitialboundaryvalueproblemiswell-posed.Inaddition,theconditionsaresuchthatthesolutionislinearlyunstable,

sowecanstudythegrowthofwaves.Weimposethat Q

˙

=

0,sothatoptionAandBareequivalent.Inthiscase,weusea standardfour-stage,fourth-orderRKmethodtointegrateintime.

We perturb the steady state by imposing a sinusoidal disturbance withwavenumber k

=

2

π

(unit 1/m)and a small amplitude.Linearstabilityanalysis[11,16] givesthefollowingangularfrequencies

ω

(unit1/s):

ω

1

=

3

.

22

+

2

.

00i

,

(60)

ω

2

=

10

.

26

−

1

.

61i

.

(61)

Thenegativeimaginarypartof

ω

2indicatestheinstabilityinthesolution.Wechoosetheinitialperturbationrelatedto

ω

1 tobezero,implyingthatasinglewavewithfrequency

ω

2results,andthelinearizedsolutionis

U

(

s

,

t

)

=

U0

+

Re





U ei(ω2t−ks)



,

(62)

where



U isobtainedbychoosingtheliquidhold-upfractionperturbationas

α

ˆ

l

=

10−3,andthencomputingthe pertur-bationsinthegasandliquidvelocityfromthedispersionanalysis[11].Theinitialconditionwhichfollowsbytakingt

=

0 s doesingeneralnotsatisfy(11) exactly,andwethereforeperformaprojectionsteptomaketheinitialconditionsconsistent, see[17],insuchawaythatthevolumetricflowratestaysexact.

Fig.3showsthesolutionatt

=

1

.

5 s computedwithN

=

40 finitevolumesand



t

=

1

/

100 s.Theamplitudehasgrown with respectto the initial condition dueto the negative imaginary part in

ω

2. At the same time, non-lineareffects are causingwavesteepening,whichwillleadtoshockformationatlatertimes.Fig.4showstheerrorinthesolutionupontime steprefinement,wheretheerroriscomputedwithrespecttoareferencesolutioncomputedwith



t

=

1

×

10−4s. Fig.5 showsthattheerrorsrelatedtothevolumeconstraint, thevolumetricflowconstraint,andtheactualvolumetricfloware satisfieduntilmachineprecision.Notethattheseerrorsaredefinedateachtimeinstantas:

volume constraint error: max

|

−_p1



_Un 1

ρg

+

Un₂

ρ

l



−

A

|,

(63)

volumetric flow constraint error: max

|

−_p1Dp



−u1



_Un 3

ρg

+

Un 4

ρ

l



|,

(64)

volumetric flow error: max

|

−_u1



_Un 3

ρg

+

Un₄

ρ

l



−

Q

(

tn

)

|,

(65)

Fig. 3. Solutionatt=1.5 s withN=40 andt=1/100 s (blue),includinginitialcondition(blackdashed).(Forinterpretationofthecoloursinthefigures, thereaderisreferredtothewebversionofthisarticle.)

Fig. 4. Fourth-order convergence of the error at t=1.5 s for the Kelvin-Helmholtz test case.

6.2. Discontinuoussolutions:rollwaves

Inthissecond testcasewe simulaterollwavesona periodicdomaininorderto (i)investigatethedifference in Q

(

t

)

betweentheproposedPF-TFMandtheoriginalTFMforthecaseofperiodicboundaryconditions,(ii)showthatthePF-TFM capturesthesame discontinuoussolutionastheTFM,and(iii)givean indicationofthe computationalspeed-upthatcan beachievedwiththePF-TFM.

Fig. 5. Errors in the volume constraint, the volumetric flow constraint and the volumetric flow remain at machine precision (t=1/100 s).

Table 2

Parametervaluesforrollwave simula-tion.

Parameter Value Unit

ρl 998 kg/m3 ρg 50 kg/m3 R 0.05 m g 9.8 m/s2 μg 1.61·10−5 Pa s μl 1·10−3 Pa s ϕ 0 deg 2·10−5 _m L 3 m

RollwavescanformwhensimulatingtheKelvin-Helmholtzinstability ofsection 6.1furtherintothenonlinearregime. The testcasethat we performhereisinspired bythe experimentalwork ofJohnson [9] andsimulations ofAkselsen[2], Holmås [7],andSanderse[15].The parametersettingsare reportedinTable 2.We firstcompute asteady statesolution: weprescribetheflowrates Qg

/

A

=

3

.

5 m

/

s, Ql

/

A

=

0

.

35 m

/

s,whichleads to

α

l

=

0

.

190 and dpds0

= −

155

.

919 Pa

/

m.The initialconditionfollowsbyapplyingasinusoidalperturbationoftheform(62),withahold-upamplitudeof0.01,k

=

2

π

/

L

and

ω

=

4

.

597

−

0

.

068i (unit1/s)chosensuchthatasingle,unstable,waveistriggered(similartotheKelvin-Helmholtztest caseinsection6.1).Oneimportantdifferencewithrespecttotheprevioustestcaseisthattheinterfacialfrictionfactor fgl intheexpressionfortheinterfacial stress(

τ

gl

=

21fgl

ρ

g

|

ug

−

ul

|(

ug

−

ul

)

) istakenas fgl

=

m

·

fg,withm

=

12

.

5.Another importantdifferenceisthatthediscretizationoftheconvectivetermsinthemomentumequationisperformedwitha first-order upwindschemeinorder toprevent numericaloscillations atthediscontinuity (shockwave),instead ofthecentral schemeusedinthefirsttestcase.

Fig.6a showsthehold-up profileatt

=

100 s computedwithboth thePF-TFM andtheTFM.At thistime instantthe wave hastravelledapproximately73timesthroughthedomain,covering adistanceofabout219m.The seeminglylarge difference inthe positionof the discontinuity (about 0.4 m) between the PF-TFM prediction and the TFM prediction is thereforein factlessthan 0.2%ofthe totaltravelled distance.Thisdifference isconsistent withtherelative difference in thevolumetricflowratepredictedby thePF-TFMandTFM,whichisalsoapproximately0.2%(seeFig.6b). Thisdifference reaches a constant value afterapproximately 50 seconds,as the rollwave then reaches a stationarysolution. Note that theTFM prediction,which doesnot requireaprescribed valueof Q

(

t

)

,should beinterpreted asbeingthemostaccurate solution.ThePF-TFMpredictioncanbeseenasahighlyaccurateapproximationtothissolutionatamuchreduced compu-tationalcost.Thedifferencebetweenthetwosolutionsdoesnotvanishupongridortime-steprefinement,sinceitiscaused by amodelassumption(namelyconstantvolumetricflow)ratherthanadiscretizationeffect.Fortunately,thisassumption isonlynecessaryforthecaseofperiodicboundaryconditions,asforothertypesofboundaryconditions,Q isknownfrom theboundaryconditionvalues.

Forthisparticularcase(N

=

320),thecomputationalcostofthePF-TFMpredictionisabout40%lowerthanoftheTFM. Forother gridsizes a similar cost reduction isobserved, asisshown inFig. 7,while thedifference inaccuracy remains negligible.Fig.7ashowsaquadraticscalingofCPUtimewithN,whichiscausedbythefactthatthesimulationsarerunat afixed CFLnumber(



t

=

1

/

N).Fig.7bshowsthattheobtainedreduction incomputationalcost isrelativelyindependent ofthenumberofgridcells.However,weshouldnotethat otherfactorscanhavean importantinfluence,suchasthetype ofpressuresolverusedintheTFM(hereadirectsolverwasused)orthetypeoffrictionmodel.Forexample,whenamore computationally expensivefrictionmodellikethe Bibergmodel[4] isused,the gainincomputational costobtainedwith thePF-TFMcomparedtotheTFMisprobablylesspronounced,becausethecontributionofthepressurePoissonsolvetothe

Fig. 6. Roll wave solutions computed with PF-TFM and TFM, with N=320 andt=1/320 s.

Fig. 7. Comparison of computational time of PF-TFM and TFM.

Table 3

ParametervaluesfortheIFPproblem. Parameter Value Unit

ρl 1003 kg/m3 ρg 1.26 kg/m3 R 0.073 m g 9.8 m/s2 μg 1.8·10−5 Pa s μl 1.516·10−3 Pa s 10−8 _m L 1000 m

totalcomputationalcostissmaller.Ontheotherhand,thecomputationsareperformedhereonasingleCPU.Muchlarger speed-upswiththePF-TFMcanbeexpectedinaparallel implementation:thePF-TFMframeworkisfullyexplicitandcan be easilyparallelized,whereas theoriginalTFMhasanimplicitcomponent(thesolutionofaPoissonequation)forwhich goodparallelscalingwillbemoredifficulttoobtain.

6.3. Hold-upwavepropagation

The finaltest caseconsiders thepropagationofahold-upwave througha horizontalpipelinecausedbyan increasing inlet gas flow rate, inspired by the test case in[14]. The parameters of the problemare described in Table 3. The ini-tial conditionsare steadystate productionwithinlet massflows ofliquid andgasof Il

=

1 kg

/

s and Ig,start

=

0

.

02 kg

/

s. Furthermore,wespecifythefollowingtime-dependentgasflowrateattheinlet,

Fig. 8. Third-order convergence of the error at t=1000 s for hold-up wave propagation. The weak BC corresponds to option A, the strong BC to option B.

Fig. 9. Hold-up wave propagation with N=40 andt=10 s. The weak BC corresponds to option A, the strong BC to option B.

whichissuchthat Q

˙

(

0

)

=

0 andfort

>

0Q

˙

(

t

)

(andhigherorderderivatives)arenonzero.Thevolumetricflows Qg

(

t

)

and

Ql aregivenby Ig

(

t

)/

ρ

g andIl

/

ρ

l,respectively.

ThesolutionswiththeweakBCprescription(optionA)andwiththestrongBCprescription(optionB)att

=

1000 s with

N

=

40 and



t

=

10 s areindistinguishable,seeFig.9a.Fig.9bshowsthat theerrorsinvolumeconstraintandvolumetric flowconstraint(equations(63)-(64))remainatmachineprecisionforbothmethods.Inaddition,optionBhastheproperty thatthenumericallycomputedvolumetricflowrateremainsexactlyequaltothespecifiedvalue.InoptionA,asmallerror ispresent,whichdecreaseswiththirdorderupontimesteprefinement.

A more quantitative assessment is obtainedby computing a referencesolution at t

=

1000 s with N

=

40,



t

=

1

×

10−2_{s, and weak boundary conditions (option A). At this small time step, the volumetric flow error stays at machine} precision alsowith weak imposition of Q and boundary conditions.Next we compute solutions at a sequence ofmuch larger time steps,



t

=

40

,

20

,

10

,

. . .

, both with weak BC (option A) and strong BC (option B). Fig. 8 shows that both optionsconvergeaccordingtothedesignorderofaccuracy(thirdorder).

7. Conclusions

Inthisarticlewehaveproposedanovelpressure-freeincompressibletwo-fluidmodelformultiphasepipelinetransport, in whichboth the pressure andthe volumeconstraint havebeen eliminated. Togetherwith a staggered-gridspatial dis-cretizationandanexplicitRunge-Kuttatime integrationmethod,thisyields anewhighly efficientmethodforsolving the incompressibletwo-fluidmodelequations.Theexplicitnaturealsoallowsforastraightforwardparallelization.Furthermore, theabsenceofthepressuremakesthemodelmoreamenabletoanalysis.Forexample,thePF-TFMoffers anewviewpoint instudyingtheissueoftheconditionalwell-posedness ofthe TFM.The ‘price’topaycomparedtothe conventionalform oftheincompressibletwo-fluid model,isthatthetime-derivativeofthe volumetricflowrateneedsto bespecifiedasan additionalinputtotheequations.Fortunately,inmostcasesforpipelinetransportdesign(e.g.steadyproduction,production ramp-up)thisvalueisknownfromtheboundaryconditions.

A natural question is whetherthe original TFM ‘pressure Poisson’ formulation [17] is still useful, given that a more efficient pressure-free alternative hasbeen proposed in thiswork. The answer isyes, fortwo reasons. First, theoriginal formulationalsoworkswhenQ is

˙

unknown(e.g.inthecaseofperiodicboundaryconditions).Second,asindicatedin[17],

theoriginal formulationhasthepotentialtobe extendedtomulti-dimensionalproblems,withapplicationtoforexample the volume-of-fluid approach. This isprobably not possible withthe proposed pressure-free approach: the fact that the pressure gradient can be expresseddirectly interms ofthe conservative variablesby using the constraintispresumably onlypossiblewhendealingwithone-dimensionalproblems.

Extensionofthepressure-freemodeltoone-dimensionalincompressibleflowwiththreeinsteadoftwo fluidphasesis straightforward. Extension withan energyequation is alsolikely to be possible,since inincompressibleflow the energy equation decouples fromthemass andmomentum equations(it can besolved oncethe massandmomentum equations havebeensolved).Ontheotherhand,extensiontocompressibleflowisnotstraightforward:incompressibleflow, differen-tiationofthevolumeconstraintleadstoapressureequationthatinvolvesthetimederivativeofthepressure.Consequently, the pressuregradient terminthe momentum equationscannot be eliminated inthesame wayasinthe incompressible model.

CRediTauthorshipcontributionstatement

B.Sanderse: Conceptualization,Methodology,Software,Writing- originaldraft. J.F.H.Buist: Writing- review &editing. R.A.W.M.Henkes: Resources,Writing- review&editing.

Declarationofcompetinginterest

The authors declare that they haveno known competingfinancial interests or personal relationships that could have appearedtoinfluencetheworkreportedinthispaper.

Appendix A. EigenvaluesofthePF-TFM

ToinvestigatetheeigenvaluesofthePF-TFM,onestudiesthetermsinthemodelthatinvolvepartialderivatives,i.e.

∂

U

∂

t

+

A

(

U

)

∂

f

(

U

)

∂

s

.

(A.1)

Wearethusinterestedintheeigenvaluesofthematrix

A

(

U

)

∂

f

(

U

)

∂

U

=

⎛

⎜

⎜

⎝

0 0 1 0 0 0 0 1

−

Alρg

(

Agg

+

u2g

)/

ρ

ˆ

−

Agρg

(

Alg

−

u_l2

)/

ρ

ˆ

2 Alρgug

/

ρ

ˆ

−

2 Agρgul

/

ρ

ˆ

Alρl

(

Agg

+

u2g

)/

ρ

ˆ

Ag

ρl

(

Alg

−

u2_l

)/

ρ

ˆ

−

2 Alρlug

/

ρ

ˆ

2 Ag

ρl

ul

/

ρ

ˆ

⎞

⎟

⎟

⎠ .

(A.2)

Theseeigenvaluesread

λ

1,2

=

Alρgug

+

Agρlul

±



AgAl

g

ρ

ˆ



ρ

−

ρg

ρl

(

u

)

2



ˆ

ρ

,

(A.3)

where



ρ

=

ρ

l

−

ρ

g and



u

=

ul

−

ug,and

λ

3,4

=

0

.

(A.4)

TheseeigenvaluesarethesameasthoseoftheTFM,seee.g.[17].

Appendix B. Pressure-freemodelfromfullydiscretetwo-fluidmodelequations

Inthisappendix wederive adiscrete PF-TFMbasedona fullydiscreteapproximation oftheTFM.Thiscorresponds to ‘optionB’.TheforwardEulerdiscretizationoftheTFMreads[17]:

Un₁₂+1

=

Un₁₂

+ 

t F_TFM,12

(

Un

),

(B.1) Un₃₄+1

=

Un₃₄

+ 

t F_TFM,34

(

Un

)

− 

t H

(

Un

)

pn

,

(B.2)



−p1Un1+1

ρg

+



−p1Un2+1

ρ

l

=

A

,

(B.3)

where H denotesthediscretizedpressuregradientoperator,withthetwocomponents

H3

(

U

)

p

= (

IpU1

/

ρg

)

 (

Gp

),

H4

(

U

)

p

= (

IpU2

/

ρl

)

 (

Gp

).

(B.4)

Ip

∈ R

Nu×Np is an interpolation matrix, G

∈ R

Nu×Np a differencing matrix(similar to Dp), and



denotes elementwise multiplication.Toarriveatthepressure-freemodelbasedonthediscretetwo-fluidmodel,wefollowexactlythesamesteps