interfaces and cavitation A numerical model for multiphase liquid–vapor–gas flows with International Journal of Multiphase Flow

ContentslistsavailableatScienceDirect

International Journal of Multiphase Flow

journalhomepage:www.elsevier.com/locate/ijmulflow

A numerical model for multiphase liquid–vapor–gas flows with interfaces and cavitation

Marica Pelanti

^a,∗

, Keh-Ming Shyue

^b

a Institute of Mechanical Sciences and Industrial Applications, UMR 9219 ENSTA ParisTech – EDF – CNRS – CEA, 828, Boulevard des Maréchaux, Palaiseau Cedex 91762, France

b Department of Mathematics and Institute of Applied Mathematical Sciences, National Taiwan University, Taipei 10617, Taiwan

a rt i c l e i nf o

Article history:

Received 2 September 2018 Revised 16 January 2019 Accepted 28 January 2019 Available online 29 January 2019 MSC:

65M08 76T10 Keywords:

Multiphase compressible flows Relaxation processes Liquid–vapor phase transition Finite volume schemes Riemann solvers

a b s t ra c t

Weareinterestedinmultiphaseflowsinvolvingtheliquidandvaporphasesofonespeciesandathird inertgaseousphase.Wedescribetheseflowsbyahyperbolicsingle-velocitymultiphaseflowmodelcom- posedofthephasicmassandtotalenergyequations,thevolumefractionequations,andthemixturemo- mentumequation.Themodelincludesstiff mechanicaland thermalrelaxationsourcetermsforallthe phases,andchemicalrelaxationtermstodescribemasstransferbetweentheliquidandvaporphasesof thespeciesthatmayundergotransition.First,wepresentananalysisofthecharacteristicwavespeeds associatedtothehierarchyofrelaxedmultiphasemodelscorrespondingtodifferentlevelsofactivation ofinfinitelyfastrelaxationprocesses,showingthatsub-characteristicconditionshold.Wethenpropose amixture-energy-consistentfinite volumemethodforthenumerical solutionofthemultiphasemodel system.Thehomogeneousportionoftheequationsissolvednumericallyviaasecond-orderwaveprop- agationschemebasedonrobustHLLC-typeRiemannsolvers.Stiff relaxationsourcetermsaretreatedby efficient numericalprocedures that exploitalgebraicequilibrium conditionsfor therelaxedstates. We presentnumericalresultsforseveralthree-phaseflowproblems,includingtwo-dimensionalsimulations ofliquid–vapor–gasflowswithinterfacesandcavitationphenomena.

© 2019ElsevierLtd.Allrightsreserved.

1. Introduction

Liquid–vaporflowsarefoundinalargevarietyofindustrialand technologicalprocessesandnaturalphenomena.Oftentheseflows involveone or more additionalinert gasphases. Forinstance, in some processes the dynamics of a liquid–vapor mixture is cou- pledtothedynamicsofdefinedregionsofathirdnon-condensable gaseouscomponent.Anexampleisgivenbyunderwaterexplosion phenomena, where a high pressure bubble of combustion gases triggers cavitation phenomena in water (Cole, 1948; Kedrinskiy, 2005).Inother casesliquid–vapormixturesmaycontainadiluted inertgascomponent, whichmayaffecttheflow features,such as infuelinjectors(Battistonietal.,2015).Weareinterestedherein thesimulationofthistype ofmultiphase flows involvingtheliq- uidandvapor phasesof onespecies andone ormore additional non-condensable gaseous phases. We describe these multiphase flowsbyahyperbolicsingle-velocitycompressibleflowmodelwith infinite-rate mechanicalrelaxation, which extends the two-phase

∗ Corresponding author.

E-mail addresses: marica.pelanti@ensta-paristech.fr (M. Pelanti), shyue@math.

ntu.edu.tw (K.-M. Shyue).

model that we have studied inprevious work Pelantiand Shyue (2014b).Thismodeliscomposedofthephasic massandtotalen- ergy equations, the volume fraction equations, and the mixture momentumequation.Themodelincludesthermalrelaxationterms toaccountforheattransferprocessesbetweenallthephases,and chemical relaxationterms to describe mass transfer betweenthe liquid and vapor phases of the species that may undergo tran- sition. Similar hyperbolic multiphase flow models withinstanta- neouspressure relaxationhave beenpreviously presented forin- stanceinPetitpasetal.(2009),Le Métayeretal.(2013),andZein etal.(2013).Afirstcontributionofourwork isa rigorousderiva- tionofthereducedpressure-relaxedmodelresultingfromthepar- ent non-equilibrium multiphase flow model with heat and mass transfertermsinthelimitofinstantaneousmechanicalrelaxation.

This is done by following the asymptotic analysistechnique em- ployed by Murrone and Guillard (2005) for the two-phase case withnoheatandmasstransfer.Moreover,we presentan original analysisofthecharacteristicwavespeedsassociatedtothehierar- chyofrelaxedmultiphasemodelscorrespondingtodifferentlevels ofactivationofinfinitelyfastmechanicalandthermo-chemicalre- laxationprocesses.Similartoresultsshownintheliteratureforthe two-phasecase(FlåttenandLund,2011;Lund, 2012;Linga,2018), https://doi.org/10.1016/j.ijmultiphaseflow.2019.01.010

0301-9322/© 2019 Elsevier Ltd. All rights reserved.

we demonstrate that sub-characteristic conditions hold, namely the speed of soundof themultiphase mixture is reducedwhen- everanadditionalequilibriumassumptionisintroduced.Then,we presenta finitevolumemethodforthe numericalsolutionofthe multiphasemodelsystembasedonaclassicalfractional steppro- cedure. The homogeneous hyperbolic portion ofthe equations is solved numericallyvia a second-order accuratewave propagation scheme, which employs a HLLC-type Riemann solver. In particu- lar,we presentherea newgeneralizedHLLC-solverbased onthe idea oftheSuliciu relaxationsolverofBouchut(2004),extending the solver that we have recently proposed in De Lorenzo et al.

(2018)forthetwo-phasecase.ThisHLLC/Suliciu-typesolverallows ustoguaranteepositivitypreservationwithasuitablechoiceofthe wave speeds. Stiff relaxationsourcetermsare treatedbyefficient numericalproceduresthatexploitalgebraicequilibriumconditions for the relaxed states, following the ideas of our work (Pelanti andShyue(2014b)).Similarapproacheshavebeenpreviouslypre- sentedintheliteratureforinstanceinLeMétayeretal.(2013).One important property of our numerical method is mixture-energy- consistencyinthesensedefinedinPelantiandShyue(2014b)),that isthemethodguaranteesconservationofthemixturetotalenergy atthediscretelevel,anditguaranteesconsistencybyconstruction ofthevaluesoftherelaxedstateswiththemixturepressurelaw.

Thispropertyisensuredthankstothetotal-energy-basedformula- tionofthemodelsystem.Wepresentseveralnumericalresultsfor three-phase flow problems, including problems involving liquid–

vapor–gas flows with interfaces and cavitation phenomena, such asunderwaterexplosiontests.

This article is organized as follows. In Section 2 we present the multiphase flow model under study. In Section 3 we derive thelimitpressure-equilibriummodelassociatedtotheconsidered multiphase flow model,andwe analyze thecharacteristic speeds oftherelaxedmodelsinthehierarchystemmingfromtheparent relaxationmodel.InSection4weillustratethenumericalmethod that we have developed to solve the multiphase flow equations.

Several one-dimensional and two-dimensional numerical experi- mentsarefinallypresentedinSection5.

2. Single-velocitymultiphasecompressibleflowmodel

We consider an inviscid compressible flow composed of N phasesthat we assume in kinematicequilibriumwith velocityu. Inthisworkwearemainlyinterestedinthree-phaseflows,N=3, nonetheless we shallpresentherea generalmultiphase flowfor- mulation. The volume fraction, density, internal energy per unit volume, and pressure of each phase will be denoted by

α

k,

ρ

k, Ek,p_k,k=1,...,N, respectively.We willdenotethetotal energy for the kth phase with E_k=Ek+

ρ

k|^u|²

2 . The saturation condition is N

k=1

α

k=1. The mixture density is

ρ

=N

k=1

α

k

ρ

k, the mix- ture internal energyis E=N

k=1

α

kEk,and themixture total en- ergy is E=N

k=1

α

kE_k=E+

ρ

^|u2^|2. Moreover, we will denote the kthspecificinternalenergywith

ε

k=Ek/

ρ

k.Mechanicalandther- maltransferprocessesareconsideredingeneralforallthephases.

We assume that one species in the mixture can undergo phase transition, so that it can exist asa vapor or a liquid phase, and mass transfer terms are accountedfor thisspecies only. We will usethesubscripts1and2todenotethe liquidandvapor phases ofthisspecies.Wedescribe theN-phaseflowunderconsideration by acompressibleflow modelthatextends thesix-equationtwo- phase flow systemthatwe studied inPelanti andShyue(2014b). The model systemis composed ofthe volume fractionequations forN− 1phases,themassandtotalenergyequationsforalltheN phases,anddmixturemomentumequations,whereddenotesthe

spatialdimension:

∂

^t

α

k+u·

∇α

k=

N j=1

Pk j, k=1,3,...,N, (1a)

∂

^t

( α

1

ρ

1

)

+

∇

·

( α

1

ρ

1u

)

=M, (1b)

∂

^t

( α

2

ρ

2

)

+

∇

·

( α

2

ρ

2u

)

=−M, (1c)

∂

^t

( α

k

ρ

k

)

+

∇

·

( α

k

ρ

ku

)

=0, k=3,...,N, (1d)

∂

^t

( ρ

^u

)

⁺

∇

^·



ρ

^uu+



_N



k=1

α

kp_k



I



=0, (1e)

∂

t

( α

1E1

)

+

∇

·

( α

1

(

^E1+p1

)

^u

)

+

ϒ

1

=−

N

j=1

p_I1jP1j+

N

j=1

Q1j+



gI+

|

^u

|

²

2

M, (1f)

∂

^t

( α

2E₂

)

⁺

∇

^·

( α

2

(

^E2+p₂

)

^u

)

⁺

ϒ

2

=−

N

j=1

pI2jP2j+

N

j=1

Q2j−



gI+

|

^u

|

²

2

M, (1g)

∂

^t

( α

kE_k

)

+

∇

^·

( α

k

(

^Ek+p_k

)

^u

)

+

ϒ

k

=−

N j=1

pIk jPk j+

N j=1

Qk j, k=3,...,N. (1h)

Thenon-conservative terms ϒk appearing inthe phasic totalen- ergyEqs.(1f)–(1h)aregivenby

ϒ

k=^→u·



Y_k

∇



_N



j=1

α

jp_j



−

∇ ( α

kp_k

)



, k=1,...,N, (1i)

whereY_k=^αk_ρ^ρk denotesthemassfractionofphase k.Inthesys- temabovePk jandQk jrepresentthevolumetransferandtheheat transfer, respectively, betweenthe phases k andj,k,j=1,...,N.

The termM indicates the mass transfer betweenthe liquid and vaporphasesindexedwith1and2.Thetransfertermsaredefined asrelaxationterms:

Pk j=

μ

k j

(

^pk− pj

)

, Qk j=

ϑ

k j

(

^Tj− Tk

)

, M=

ν (

^g2− g1

)

, (2) whereT_k denotes thephasic temperature,g_kthe phasicchemical potential,andwherewehaveintroducedthemechanical,thermal andchemical relaxation parameters

μ

k j=

μ

jk≥ 0,

ϑ

k j=

ϑ

jk≥ 0, and

ν

=

ν

12=

ν

21≥ 0, respectively. Note that: Pkk=0, Qkk=0, Pk j=−PjkandQk j=−Qjk.Thequantitiesp_{Ik j}=p_Ijk areinterface pressures and g_I is an interface chemical potential. We shall as- sume that all mechanical relaxation processes are infinitely fast,

μ

k j=

μ

jk≡

μ

→+∞,so that mechanicalequilibriumis attained instantaneouslybetweenallthephases.Indeedhere,followingthe same idea of Saurel et al. (2008, 2009) and Pelanti and Shyue (2014b)),theparent non-equilibriummultiphaseflowmodelwith instantaneouspressurerelaxationisusedtoapproximatesolutions tothelimitingpressure-equilibriumflowmodel,whichisthephys- icalflow model ofinterest. Concerningthermal andchemical re- laxation, following the simple approach of Saurel et al. (2008), we consider in this work that these processes are either inac- tive,

ϑ

k j=0,

ν

=0,ortheyactinfinitelyfast,

ϑ

k j→+∞,

ν

→+∞. Heatandmasstransfermaybeactivatedatselectedlocations,for

instanceatinterfacesforaphase pair(k,j),identifiedbymin(

α

k,

α

j)>

^,where

^{isatolerance.}

Theclosureofthesystem(1)isobtainedthroughthespecifica- tionofan equation ofstate (EOS)foreach phase p_k=p_k(Ek,

ρ

k), T_k=T_k(^pk,

ρ

k)^{. For the numerical model here we will adopt the} widely used stiffened gas (SG) equation of state (Menikoff and Plohr,1989):

p_k

(

^Ek,

ρ

k

)

⁼

( γ

k− 1

)

^Ek−

γ

k



k−

( γ

k− 1

) η

k

ρ

k, (3a)

T_k

(

^pk,

ρ

k

)

^{= pk+}



k

κ

vk

ρ

k

( γ

k− 1

)

^{, (3b)}

where

γ

k,ϖk,

η

kand

κ

vk areconstantmaterial-dependentparam- eters.Inparticular,

κ

vkrepresentsthespecificheatatconstantvol- ume.Thecorrespondingexpressionforthephasicentropyis s_k=

κ

vklog

(

^T_k^γk

(

^pk+



k

)

^−(γk−1)

)

+

η

k^, (3c) where

η

^k=constant,andg_k=h_k− Tks_k,withh_kdenotingthepha- sicspecific enthalpy. The parameters for the SGEOS for theliq- uidandvapor phasesofthe speciesthat may undergotransition are determined so that the theoretical saturation curves defined by g₁=g₂ fit the experimental onesfor the considered material (LeMétayeretal.,2004). Themixturepressurelawforthemodel withinstantaneous pressurerelaxationisdetermined bythemix- tureenergyrelation

E=

N

k=1

α

kEk

(

p,

ρ

k

)

, (4)

wherewehaveusedthemechanicalequilibriumconditionsp_k=p, forallk=1,...,N,inthephasicenergylawsEk(^pk,

ρ

k)^.Notethat fortheparticularcaseoftheSGEOS,anexplicitexpressionofthe mixturepressurecanbeobtainedfrom(4).

Sinceherewewillconsiderrelaxationparameterseither=0or

→+∞,a specificationoftheexpression fortheinterface quanti- tiesp_Ikj,g_Iisnotneeded.Nevertheless,letusremarkthatthedef- initionof these interface quantities must be consistent withthe secondlawofthermodynamics,whichrequiresanon-negativeen- tropyproductionforthemixture.Theequationforthemixtureto- talentropyS=

ρ

^s,s=N

k=1Y_ks_k,isfoundas:

∂

tS+

∇

·

(

Su

)

=HP+HQ+HM, (5a) where

HP =

N

k=1

N j=1

pk− pIk j

T_k Pk j, HQ=

N

k=1

N j=1

1 T_kQk j, HM=

_g

I− g1

T1 −gI− g2

T2



M. (5b)

For consistency of the multiphase model (1)with the second lawofthermodynamicsweneedHP+HQ+HM≥ 0.Byfollowing theargumentsinFlåttenandLund(2011),onecaninferthefollow- ingsufficientconsistencyconditionsontheinterfacequantities:

p_{Ik j}∈[min

(

^pk,p_j

)

,max

(

^pk,p_j

)

^]

and gI∈[min

(

^g1,g2

)

^,max

(

^g1,g2

)

^{]. (6)}

Themodel(1)ishyperbolic andtheassociatedspeed ofsoundc_f (non-equilibriumorfrozensoundspeed)isdefinedby

c²_f =

 ∂

^pm

∂ρ

sk,Yk,αk, k=1,...,N

, (7)

wherewehaveintroducedthemixturepressure pm

( ρ

,s1,...,sN,Y1,...,YN−1,

α

1,...,

α

N−1

)

=

N

k=1

α

kpk

sk,

ρ α

^Ykk



.

(8) Fromthisdefinition,bynoticingthat

 ∂

^pk

∂ρ

sk,Yk,αk

=

 ∂

^pk

∂ρ

k

sk,Yk,αk

 ∂ρ

k

∂ρ

sk,Yk,αk

=c²_kY_k

α

k

, (9)

weobtaintheexpression:

c_f=



N k=1

Y_kc²_k, (10)

wherec_k=

(^∂_∂ρ^pk

k)s_k isthespeedofsoundofthephasek,which can be expressed as c_k=



kh_k+

χ

k, where k=(

∂

^pk/

∂

Ek)ρ_k (Grüneisencoefficient),and

χ

k=(

∂

^pk/

∂ ρ

k)Ek.

3. Hierarchyofmultiphaserelaxedmodelsandspeedofsound

Byconsideringdifferentlevelsofactivationofinstantaneousre- laxationprocesseswecanestablishfromthemodel(1)ahierarchy of hyperbolic multiphase flow models. Here in particularwe de- rivethe expressionofthe speedofsoundforthe relaxedmodels in this hierarchy, similar to Flåtten and Lund (2011) and Flåtten etal.(2010).

3.1. p-Relaxedmodel

Intheconsidered limit ofinstantaneousmechanicalrelaxation

μ

k j≡

μ

→ +∞, the model system (1) reduces to a hyperbolic single-velocitysingle-pressuremodel,whichisa generalizationof the five-equation two-phase flow model of Kapila et al. (2001). Thereducedpressureequilibriummodel,whichweshallalsocall p-relaxedmodel,canbe derived bymeans ofasymptoticanalysis techniques,cf.inparticularMurroneandGuillard(2005).Weshow thederivationfortheone-dimensionalcaseinAppendixA.Denot- ing withp the equilibriumpressure, we obtain thefollowing re- laxedsystem,composedof2N+d equations:

∂

^t

α

1+u·

∇α

1=K₁

∇

^{· u+}

ρ ^

^1c121

N

j=2

Q1j

−

α

1

ρ

^c2p

ρ

1c²₁

N

j,i=1 i> j

Qji

 _

j

ρ

jc²_j −



i

ρ

ic_i²

+

ρ

^c2p

ρ

1c²₁

×



(

1

(

^gI− h1

)

+ c²₁

)

^N

j=2

α

j

ρ

jc²_j +

(

2

(

^gI− h2

)

+ c²₂

) α

¹

ρ

2c²₂



M, (11a)

∂

^t

α

k+u·

∇α

k=K_k

∇

^{· u+}

ρ ^

kc^k2_k

N

j=1 j=k

Qk j

−

α

k

ρ

^c2p

ρ

kc²_k

N

j,i=1 i> j

Qji

 _

j

ρ

jc²_j −



i

ρ

ic_i²

+

ρ

^c2p

α

k

ρ

kc²_k

×

 

2

(

^gI− h2

)

+c²₂

ρ

^2c2²

−



1

(

^gI− h1

)

+c²₁

ρ

^1c21

M, k=3,...,N, (11b)

∂

t

( α

1

ρ

1

)

+

∇

·

( α

1

ρ

1u

)

=M, (11c)

∂

^t

( α

2

ρ

2

)

+

∇

·

( α

2

ρ

2u

)

=−M, (11d)

∂

^t

( α

k

ρ

k

)

⁺

∇

^·

( α

k

ρ

ku

)

^{=0, k=3,...,N, (11e)}

∂

^t

( ρ

^u

)

+

∇

^·

( ρ

^uu+pI

)

=0, (11f)

∂

^tE+

∇

^·

((

^E+p

)

^u

)

=0, (11g) where

Kk=

ρ

^c2p

α

k

N

j=1 j=k

α

j



1

ρ

kc²_k − 1

ρ

jc²_j

=

α

k

 ρ

^c2p

ρ

kc²_k− 1

. (12)

Intherelationsabovewehaveintroducedthepressureequilibrium speed of soundcp (ageneralization of Wood’ssound speed),de- finedby

c²_p=

 ∂

^p

∂ρ

s1,...,sN,Y1,...,YN

, (13)

fromwhichweobtaintheexpression:

cp=

 ρ

^N

k=1

α

k

ρ

kc²_k



−¹2

. (14)

As we mentioned above, the pressure equilibrium model (1) is indeed the physical flow model of interest. Similar to the two- phasecase(Saureletal.,2009;Zeinetal.,2010;PelantiandShyue (2014b)), thenon-equilibriummodel (11)withinstantaneous me- chanicalrelaxationisconvenienttoapproximatenumericallysolu- tionstothep-relaxedmodel.

Remark1. Forthetwo-phasecaseN=2thep-relaxedmodel(11) hasaformanalogoustothepressure-equilibriummodelpresented bySaurel etal.(2008),nonethelesswe remarkadifference inthe expressionofmasstransfertermappearinginthevolumefraction equation.Theequationfor

α

1obtainedfrom(11)forN=2canbe writtenas:

∂

^t

α

^1+u·

∇α

^1=K1

∇

^{· u+}

ζ

 

1

α

^{1 +}

^ α

²²

Q +

ζ

 

1

(

^gI− h1

)

+c²₁

α

1

+



2

(

^gI− h2

)

+c²₂

α

2

M, (15)

where K₁=

ζ

(

ρ

2c²₂−

ρ

1c²₁) ^and

ζ

=_α ^α1α2

2ρ1c₁²+α1ρ2c²₂. The equation for the volume fraction

α

1 of the relaxed pressure-equilibrium modelreportedinSaureletal.(2008)is:

∂

^t

α

1+u·

∇α

1=K₁

∇

· u+

ζ

 

1

α

1 +



2

α

2

Q+

ζ



c²₁

α

1 +c²₂

α

2

M. (16) Weobservethatthetwoformulationsareequivalentonlywiththe followingdefinitionoftheinterfacechemicalpotentialg_I: g_I=

α

2



1h1+

α

1



2h2

α

2



1+

α

1



2 . (17)

Thisdefinitioningeneraldoesnot satisfythe sufficientcondition forentropyconsistency (6).Nevertheless,let usnote thatthenu- merical model in Saurel et al. (2008) considers either no mass transferorinfinite-ratemasstransfer, sothat thefactormultiply- ingMin(16)hasnoinfluenceinthesespecificcircumstances.

Remark 2. In our previous work (Pelanti and Shyue (2014b)) an additionaltermoftheformM/

ρ

_I^{waswritteninthevolumefrac-}

tion equation of the six-equation two-phase flow model corre- spondingto(1)forN=2,with

ρ

_I representinganinterface den- sity.SimilartoFlåttenandLund(2011),thistermisnotincludedin thepresentmultiphasemodel(1).The purposeofthetermM/

ρ

_I

in Pelanti and Shyue (2014b)) was to indicate the influence of themasstransferprocessontheevolutionofthevolumefraction.

Nonetheless,therigorousderivationofthepressure-relaxedmodel (11) fromthesystem(1) revealsthat indeedmasstransfer terms affect

α

kvia thepressure relaxationprocess,aswe observefrom thecontribution ofM appearingin (11a)and(11b)(and (15)for the caseN=2). Note that neglecting the term M/

ρ

I inthe six- equationtwo-phase modelofPelantiandShyue(2014b) doesnot affect the numerical model and the numerical results presented there,since

ν

= 0 or

ν

→+∞, andthe numericalprocedure for treatinginstantaneouschemicalrelaxationconsistsinimposingdi- rectlyalgebraicthermodynamicequilibriumconditions.

3.2.pT-relaxedmodels

Assuming instantaneous mechanical equilibrium

μ

jk≡

μ

→ +∞ for all the phases and thermal equilibrium

ϑ

k j≡

ϑ

→ +∞ for M phases, 2≤ M≤ N, we obtain a hyperbolic relaxed system of2N− M+ 1 + d equations characterized by the speedof sound c_pT,M,definedby

1 c_pT,M² =

 ∂

^p

∂ρ

M

k=1Yksk,sM+1,...,sN,Y1,...,YN

, (18)

Fromthisdefinitionweobtaintheexpression:

1 cpT,M2 = 1

cp2+

ρ

^T

M k=1C_pk

M−1

k=1

C_pk

M

j=k+1

Cp j

 

j

ρ

jc²_j −



k

ρ

kc²_k

2

, (19)

where T denotes the equilibrium temperature, C_pk=

α

k

ρ

k

κ

pk,

κ

pk=(

∂

^hk/

∂

^Tk)^p_k ^{(specificheat atconstant pressure),andwe re-} call k=(

∂

^pk/

∂

Ek)ρk. Let us note that in the particular case of thermalequilibriumforallthephases,M=N,thereducedsingle- pressuresingle-temperaturepT-relaxedmultiphasemodel hasthe conservativeform:

∂

t

( α

1

ρ

1

)

+

∇

·

( α

1

ρ

1u

)

=M, (20)

∂

^t

( α

²

ρ

²

)

⁺

∇

^·

( α

²

ρ

^2u

)

^{=−M, (21)}

∂

t

( α

k

ρ

k

)

+

∇

·

( α

k

ρ

ku

)

=0, k=3,...,N, (22)

∂

^t

( ρ

^u

)

⁺

∇

^·

( ρ

^uu+pI

)

^{=0, (23)}

∂

tE+

∇

·

((

^E+p

)

^u

)

=0. (24) The two-phase (N=2) version of this modelwas considered for instanceinLundandAursand(2012),andmorerecentlyinSaurel etal.(2016).

3.3.pTg-relaxedmodels

We now assume instantaneous mechanical equilibrium

μ

jk≡

μ

→+∞ for all the phases, thermal equilibrium

ϑ

k j≡

ϑ

→+∞ for M phases, 2≤ M≤ N, and, additionally, we consider instanta- neouschemical relaxationbetweentheliquidandvaporphases1 and2,

ν

→+∞.Weconsiderthat atleasttheliquid–vaporphase pairisinthermalequilibrium.Withthesehypothesesweobtaina hyperbolicrelaxedsystemof2(^N− M+ 1)+ d equationscharacter- izedbyaspeedofsoundc_pTg,M,definedby

cpTg,M2=

 ∂

^p

∂ρ

M

k=1Yksk,sM+1,...,sN,Y3,...,YN

, (25)

fromwhichweobtain

1

cpTg,M2 = 1

cpT,M2 +

ρ

^T

M k=1C_pk



M k=1



kC_pk

ρ

kc²_k −1 T



dT dp

sat

M

k=1

C_pk



2

, (26)

wherewehaveintroduced thederivatives(dT/dp)sat evaluated on theliquid–vaporsaturationcurve.Asexpected(cf.e.g.Stewartand Wendroff,1984), analogouslyto the two-phase case(Flåtten and Lund, 2011), it is easy to observe from (14), (19), and (26) that sub-characteristic conditionshold, namelythe speed of sound of the N-phase mixture is reduced whenever an additional equilib- riumassumptionisintroduced:

cpTg≡ cpTg,N≤ cpTg,M, cpT ≡ cpT,N≤ cpT,M, and

cpTg<cpT <cp<c_f. (27)

Letusnotethat in theparticularcaseof thermalequilibriumfor allthe phases,M=N,the reducedpTg-relaxedmultiphase model corresponds tothe well known Homogeneous EquilibriumModel (HEM)(Stewart and Wendroff,1984), composed ofthe conserva- tionlawsforthemixture density

ρ

^{, themixturemomentum}

ρ

^u, andthe mixturetotal energy E. The derivation of theexpression ofthespeedofsoundforthe consideredhierarchyofmultiphase flowmodelsisdetailedinAppendixB.Weconcludethissectionby showinginFig.1thebehaviorofthesoundspeedfordifferentlev- elsofactivation ofinstantaneous mechanical,thermalandchemi- calrelaxationforathree-phasemixturemadeofliquidwater,wa- tervaporandair(non-condensablegas).Hereweplotthespeedof soundversus thevolumefractionofthetotalgaseous component

α

^gv=

α

^v+

α

^{g fora fixedratio}

α

^g/

α

^v=0.5,where here

α

^{v is the}

vaporvolumefraction,and

α

gisthenon-condensablegasvolume fraction.The referencepressure is p=10⁵ Pa, and the reference temperatureisthecorresponding saturationtemperature. Thepa- rametersusedfortheequationsofstateofthephasesarethesame asthoseofthecavitation tubeexperimentinSection5.2(Experi- ment5.2.1).

4. Numericalmethod

We focus now on the numerical approximation of the multi- phasesystem (1),whichwe can write in compactvectorial form denotingwithq∈R^3N−1+d thevectoroftheunknowns:

∂

^tq+

∇

· F

(

^q

)

+

ς (

^q,

∇

^q

)

=

ψ

μ

(

^q

)

+

ψ

ϑ

(

^q

)

+

ψ

ν

(

^q

)

, (28a)

Fig. 1. Speed of sound for a three-phase mixture made of liquid water, water vapor and air versus the total gaseous volume fraction αvg = αv + αg . We use the sub- scripts l,v,g to indicate the liquid phase, the vapor phase and the non-condensable gas phase, respectively. c f , c ^pT = c pT,3 , c pT(jk) = c pT,2 , c pTg = c pTg,3 , c pT(jk)g = c pTg,2 are the speeds defined in (10), (14), (19) , and (26) . Here the notation T ( jk ), j, k = l , v , g , specifies the two phases for which thermal equilibrium is assumed (for instance c pT(lv) denotes the speed of sound for a mixture characterized by pressure equilib- rium for all the phases and thermal equilibrium for the liquid and vapor pair only).

q=

⎡

⎢ ⎢

⎢ ⎢

⎢ ⎢

⎢ ⎢

⎢ ⎢

⎢ ⎢

⎢ ⎢

⎢ ⎢

⎢ ⎢

⎢ ⎢

⎣ α

1

α

³

..

α

.N

α

1

ρ

1

α

²

ρ

²

..

α

^N.

ρ

^N

ρ

^u

α

1E₁

α

^2E2

..

α

^N.^EN

⎤

⎥ ⎥

⎥ ⎥

⎥ ⎥

⎥ ⎥

⎥ ⎥

⎥ ⎥

⎥ ⎥

⎥ ⎥

⎥ ⎥

⎥ ⎥

⎦

, F

(

^q

)

=

⎡

⎢ ⎢

⎢ ⎢

⎢ ⎢

⎢ ⎢

⎢ ⎢

⎢ ⎢

⎢ ⎢

⎢ ⎢

⎢ ⎢

⎢ ⎢

⎣

0 0 .. .

α

¹0

ρ

^1u

α

2

ρ

2u

..

α

N

ρ

.Nu

ρ

^uu+



_N

k=1

α

kp_k

 α

1

(

^E1+p1

)

^u I

α

2

(

^E2+p₂

)

^u ..

α

^N

(

^EN.^+pN

)

^u

⎤

⎥ ⎥

⎥ ⎥

⎥ ⎥

⎥ ⎥

⎥ ⎥

⎥ ⎥

⎥ ⎥

⎥ ⎥

⎥ ⎥

⎥ ⎥

⎦

,

ς (

^q,

∇

^q

)

⁼

⎡

⎢ ⎢

⎢ ⎢

⎢ ⎢

⎢ ⎢

⎢ ⎢

⎢ ⎢

⎢ ⎢

⎢ ⎢

⎢ ⎢

⎢ ⎢

⎣

u·

∇α

1

u·

∇α

³

.. .

u·

∇α

^N

0 0 .. . 0

ϒ

01

ϒ

2

..

ϒ

.N

⎤

⎥ ⎥

⎥ ⎥

⎥ ⎥

⎥ ⎥

⎥ ⎥

⎥ ⎥

⎥ ⎥

⎥ ⎥

⎥ ⎥

⎥ ⎥

⎦

, (28b)

ψ

μ

(

^q

)

⁼

⎡

⎢ ⎢

⎢ ⎢

⎢ ⎢

⎢ ⎢

⎢ ⎢

⎢ ⎢

⎢ ⎢

⎢ ⎢

⎢ ⎢

⎢ ⎢

⎢ ⎣

N j=1P1j

N j=1P3j

..

_N.

j=1PN j

0 0 .. . 0 0

−_N

j=1pI1jP1j

−_N

j=1pI2jP2j

.. .

−N

j=1pIN jPN j

⎤

⎥ ⎥

⎥ ⎥

⎥ ⎥

⎥ ⎥

⎥ ⎥

⎥ ⎥

⎥ ⎥

⎥ ⎥

⎥ ⎥

⎥ ⎥

⎥ ⎦

,

ψ

θ

(

^q

)

⁼

⎡

⎢ ⎢

⎢ ⎢

⎢ ⎢

⎢ ⎢

⎢ ⎢

⎢ ⎢

⎢ ⎢

⎢ ⎢

⎢ ⎢

⎢ ⎢

⎢ ⎣

0 0 .. . 0 0 0 .. . 0

N0

j=1Q1j

N j=1Q2j

..

N .

j=1QN j

⎤

⎥ ⎥

⎥ ⎥

⎥ ⎥

⎥ ⎥

⎥ ⎥

⎥ ⎥

⎥ ⎥

⎥ ⎥

⎥ ⎥

⎥ ⎥

⎥ ⎦

,

ψ

ν

(

^q

)

⁼

⎡

⎢ ⎢

⎢ ⎢

⎢ ⎢

⎢ ⎢

⎢ ⎢

⎢ ⎢

⎢ ⎢

⎢ ⎢

⎢ ⎢

⎢ ⎢

⎢ ⎢

⎣

0 0 .. . 0

−MM .. . 0

0

g_I+^|u₂^|2



M

−

gI+^|u₂^|2



M .. . 0

⎤

⎥ ⎥

⎥ ⎥

⎥ ⎥

⎥ ⎥

⎥ ⎥

⎥ ⎥

⎥ ⎥

⎥ ⎥

⎥ ⎥

⎥ ⎥

⎥ ⎥

⎦

, (28c)

with ϒk(q,

∇

^{q) defined in (1i). Above we have put into evi-}

dencethe conservative portionof thespatialderivative contribu- tions in thesystemas

∇

· F(^q),andwe have indicated thenon- conservative term as

ς

^(q,

∇

^{q). The source terms}

ψ

μ(q),

ψ

θ(q),

ψ

ν(q)containmechanical,thermalandchemicalrelaxationterms, respectively,asexpressedin(2).

To numerically solve the system (28) we use the same tech- niques that we have developed for the two-phase model in Pelanti andShyue (2014b).Afractional step methodisemployed, where we alternate between the solution of the homogeneous system

∂

^tq+

∇

· F(^q)+

ς

(^q,

∇

^q)=0 and the solution of a se- quence of systems of ordinary differential equations (ODEs) that take into account the relaxationsource terms

ψ

μ,

ψ

ϑ, and

ψ

ν. AsinPelantiandShyue(2014b),theresultingmethodismixture- energy-consistent, inthesense that (i) itguaranteesconservation atthediscretelevelofthemixturetotalenergy;(ii)itguarantees consistencybyconstructionofthevaluesoftherelaxedstateswith the mixturepressure law.The method hasbeenimplemented by usingthelibrariesoftheclawpacksoftware(LeVeque).

4.1. Solutionofthehomogeneoussystem

To solvethe hyperbolic homogeneous portion of(28) we em- ploy the wave-propagation algorithms of LeVeque (2002, 1997), which are aclass of Godunov-typefinite volumemethods to ap- proximate hyperbolicsystemsofpartialdifferentialequations.We shallconsiderhereforsimplicitytheone-dimensional caseinthe xdirection(d=1),andwereferthereadertoLeVeque(2002)for a comprehensivepresentationofthesenumericalschemes. Hence we consider here the solution of the one dimensional system

∂

^tq+

∂

^xf(^q)+

ς

(^q,

∂

^xq)=0,q∈R^3N (asobtainedbysettingu=u and

∇

=

∂

^{x in(28)).We assumea gridwithcellsofuniformsize}

^{x,andwedenotewithQ}iⁿ theapproximatesolutionof thesys- tematthe ith celland attime tⁿ, i∈Z, n∈N.The second-order wavepropagationalgorithmhastheform

Q_iⁿ⁺¹=Q_iⁿ−



^t



^x

(

^A+



^Qi−1/2+A⁻



^Qi+1/2

)

⁻



^t



^x

(

^Fi^h+1/2− F_i^h_−1/2

)

^. (29) Here A^∓^Qi+1/2 are the so-called fluctuations arising from Rie- mannproblemsatcellinterfaces (ⁱ+1/2)^{betweenadjacentcells} iand (ⁱ+1),andF_i^h_+1/2 arecorrection termsfor(formal)second- order accuracy. To define the fluctuations, a Riemann solver (cf.GodlewskiandRaviart,1996;Toro,1997;LeVeque,2002)must beprovided.Thesolutionstructuredefinedbyagivensolverfora Riemann problem withleft and rightdata q_ and qr can be ex- pressed in general by a set of M waves W^l and corresponding speedss^l,M3N.Forexample,fortheHLLC-typesolverdescribed belowM=3.Thesumofthe wavesmustbe equaltothe initial jumpinthevectorqofthesystemvariables:



^q≡ qr− q_=^M

l=1

W^l. (30)

Moreover,foranyvariableofthemodelsystemgovernedbyacon- servativeequation theinitialjump intheassociatedfluxfunction mustbe recoveredby the sumofwavesmultipliedby thecorre- sponding speeds. In the considered modelthe conserved quanti- tiesare

α

k

ρ

k,k=1,...,N,and

ρ

^{u,thereforeinordertoguarantee}

conservationweneed:



^{f(ξ )}≡ f^{(ξ )}

(

^qr

)

− f^{(ξ )}

(

^q

)

=^M

l=1

s^lW^{l(ξ )} (31)

for

ξ

=N,. . .,2N, where f^(ξ) is the

ξ

^{th component of the flux}

vector f, andW^{l(ξ )} denotes the

ξ

^{th component ofthe lth wave,}

l=1,. . .,M. It is clear that conservation of the partial densities ensuresconservationofthemixturedensity

ρ

=N

k=1

α

k

ρ

k.Inad- dition,we must ensure conservationof themixture total energy,



^fE≡ fE

(

^qr

)

− fE

(

^q

)

=

M l=1

s^l

N

k=1

W^l(2N+k), (32) where f_E=u(^E+N

k=1

α

kp_k) ^{is the flux function associated to} the mixture total energy E. Once the Riemann solution struc- ture

{

W_i^l_+1/2,s^l_i+1/2

}

l=1,...,M arising at each celledge x_i+1/2 isde- fined through a Riemann solver, the fluctuations A^∓^Qi+1/2 and thehigher-order(second-order)correctionfluxesF_i^h_+1/2in(29)are computedas

A^±



^Qi+1/2=

M l=1

(

^sli+1/2

)

^±W_i^l_+1/2, (33)

wherewe haveused the notations⁺=max(^s,0),s⁻=min(^s,0), and

F_i^h_+1/2=1 2

M l=1



s^l_i+1/2

 

1−



^t



^x



s^l_i+1/2



W_i^l₊₁^h_/2, (34)

whereW_i^l₊₁^h_/2 areamodifiedversionofW_i^l_+1/2obtainedbyapply- ingtoW_i^l_+1/2alimiterfunction(cf.LeVeque,2002).

One difficulty in the solution of the homogeneous portion of themultiphasesystem(28)isthepresenceofthenon-conservative products ϒk in the phasic energy equations. Although a discus- sion ofthe treatment of non-conservative terms is not the main focusofthepresentwork,it isimportanttorecall theassociated issuesandchallenges.Itiswellknownthatafirstdifficultyofnon- conservative hyperbolic systems is the lack of a notion of weak