compressible multiphase flows with moving interfaces High fidelity discontinuity-resolving reconstruction for Journal of Computational Physics

Contents lists available atScienceDirect

Journal of Computational Physics

www.elsevier.com/locate/jcp

High fidelity discontinuity-resolving reconstruction for compressible multiphase flows with moving interfaces

Xi Deng

^a

, Satoshi Inaba

^a

, Bin Xie

^a

, Keh-Ming Shyue

^b

, Feng Xiao

^a,∗

aDepartmentofMechanicalEngineering,TokyoInstituteofTechnology,2-12-1Ookayama,Meguro-ku,Tokyo,152-8550,Japan bDepartmentofMathematics,NationalTaiwanUniversity,Taipei106,Taiwan

a r t i c l e i n f o a b s t r a c t

Articlehistory:

Received26October2017

Receivedinrevisedform23February2018 Accepted25March2018

Availableonline28March2018

Keywords:

Compressiblemultiphaseflows Five-equationmodel Interfacecapturing THINCreconstruction BVDalgorithm Finitevolumemethod

We present in this work a new reconstruction scheme, so-called MUSCL-THINC-BVD scheme, to solve the five-equation model for interfacial two phase flows. This scheme employsthetraditionalshockcapturingMUSCL(MonotoneUpstream-centeredSchemesfor ConservationLaw)schemeaswellastheinterfacesharpeningTHINC(TangentofHyperbola for INterface Capturing)schemeas two building-blocksof spatialreconstruction onthe BVD(boundaryvariationdiminishing) principlethatminimizes thevariations(jumps)of thereconstructedvariablesatcellboundaries,andthuseffectivelyreducesthedissipation errorinnumericalsolutions.TheMUSCL-THINC-BVDschemeisimplementedtothevolume fractionandotherstatevariablesunderthesamefinitevolumeframework,whichrealizes theconsistencyamongvolumefractionandotherphysicalvariables.Numericalresultsof benchmarktestsshowthatthepresentmethodisabletocapturethematerialinterfaceas awell-definedsharpjumpinvolumefraction,andobtainnumericalsolutionsofsuperior quality incomparison tootherexistingmethods. Theproposed schemeisasimple and effectivemethodofpracticalsignificanceforsimulatingcompressibleinterfacialmultiphase flows.

©^{2018ElsevierInc.Allrightsreserved.}

1. Introduction

Compressible multiphase fluid dynamics is one of active and challenging research areas ofgreat importance in both theoretical studies and industrial applications.For example, shock/interface interactions are thought to be crucial to the instability andevolution ofmaterialinterfaces that separate differentfluidsasobserved ina widespectrum ofphenom- ena [2].The material interfacesgreatly complicate thephysics andmake problemsformidably difficult foranalytical and experimentalapproachesinmanycases,wherenumericalsimulationturnsouttobethemosteffectiveapproachtoprovide quantitativeinformationtoelucidatethefundamentalmechanismsbehindthecomplexphenomenaofmultiphaseflows.

Incomparison tothe computationof single phaseflow, development ofnumerical methodsformultiphase flow faces morechallengingtasks.Themajorcomplexitycomesfromthemovinginterfacesbetweendifferentfluidsthatusuallyasso- ciatewithstrongdiscontinuities,singularforcesandphasechangesinsomecases.Giventhenumericalmethodsdeveloped formultiphaseincompressibleflowswithinterfaceshavingbeenreaching arelativelymaturestage, thenumericalsolvers forcompressibleinterfacialmultiphaseflowsare apparentlyinsufficient.Forincompressiblemultiphaseflowswithmoving

*

Correspondingauthor.

E-mailaddress:xiao.f.aa@m.titech.ac.jp(F. Xiao).

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

0021-9991/©^{2018ElsevierInc.Allrightsreserved.}

interfaces where the density andother physicalproperties, e.g. viscosity and thermalconductivity, are constant ineach fluid,theone-fluidmodel[1] canbeimplementedinastraightforwardmannerwithanassumptionthatthephysicalfields change monotonicallyacrosstheinterface region.Thus, providedwithan indicationfunction whichidentifies themoving interface, one can uniquely determinethe physicalproperty fieldsfor thewhole computational domain.Some indication functions, suchasvolumeoffluid (VOF)function[3–5] andlevelsetfunction[6–8],havebeenproposedandprovedtobe abletocapturemovinginterfaceswithcompactthicknessandgeometricalfaithfulnessifsolvedby advancednumericalal- gorithms.However,substantialbarrierexistswhenimplementingtheone-fluidmodeltocompressibleinterfacialmultiphase flows.

Thenewdifficultieswefacewhenapplyingtheone-fluidmodel¹ tocompressibleinterfacialmultiphaseflowslieintwo aspects:

(I) Densityandenergyincompressibleflowhavetobesolvedalongwiththeindicationfunction,andspecialformulations arerequiredtomaintainthephysicalconsistencywhichresultsinabalancedstateamongallvariablesfortheinterface cellwhereawell-definedinterfacefallsin;

(II) Thenumericaldissipationintheso-calledhigh-resolutionschemesdesignedforsolvingsinglephasecompressibleflows involvingshockwavestendstosmearoutdiscontinuitiesincludingthematerialinterfacesinnumericalsolutions,which isfataltosimulationsofinterfacialmultiphaseflowseveniftheschemescanproduceacceptableresultsinsinglephase cases.

Forissue(I)mentioned above,mixingoraveragingmodelsthat consistofEulerorNavier–Stokesequationsalongwith the equations ofinterface-indication functions have been derived andwidely used asan efficient approximation to the state of the interface cell where two or more species co-exist. A simple single-fluid model was reported in [9–11] for interfacial multiphasecompressibleflowsusingeitherexplicittime marchingorsemi-implicitpressure-projection solution procedure.Thelatterresultsinaunifiedformulationforsolvingbothcompressibleandincompressiblemultiphaseflows.As theprimitivevariablesare solvedinthesemodels,theconservationpropertiesarenotguaranteed,andthusmightnotbe suitableforhigh-Machflowsinvolvingshockwaves.Conservativeformulations,whichhavebeenwell-establishedforsingle phase compressibleflows withshockwaves,howevermayleadtospurious oscillationsinpressureorotherthermalfields [12,13]. Itwas foundthat specialtreatments arerequiredintransportingthematerial interfaceandmixing/averagingthe statevariablestofindthemixedstateoffluidsintheinterfacialcellthatsatisfiespressurebalanceacrossmaterialinterface formultiplepolytropicandstiffgases[14–18],vanderWaals[19] andMie–Grüneisenequationsofstate(EOS)[20].Amore general five-equationmodel [21] was developedfora wide rangefluids. Thesemodels apply tomultiphase compressible flows witheitherspreadorsharpinterfaces.See[22] forarecentreviewonthe modelsofthissort. Wemakeuseofthe five-equation model inthe presentwork asthe PDE(partial differential equation) setto develop ournumerical method, whichcanbeappliedtootherextendedsystemaswell.

Provided the SEF modelswith some desiredproperties,such ashyperbolicity, conservationandwell-balanced mixing closure without spurious oscillations in thermal variables, we can in principle implement numerical methods for single phase compressible flow (e.g.the standardshock-capturing schemes)to solvethesemultiphase models. TVD(Total Vari- ation Diminishing) schemes, such as the MUSCL (Monotone Upstream-centered Schemes for Conservation Law) scheme [23], can solve discontinuities without numerical oscillations, which is of paramount importance to ensure the physical fields to be bounded and monotonic in the transition region. However, TVD schemes suffer from excessive numerical dissipation, which brings the problem (II) listed above to us. The intrinsic numerical dissipation smears out the flow structures including the discontinuities in mass fraction or volume fraction which is used to represent the material in- terfaces. Consequently, material interfaces are continuously blurred and smeared out, which is not acceptable in many applications, especially for the simulations that need long-termcomputation. As a remedy, using higher order schemes, like WENO(Weighted Essentially Non-Oscillatory) scheme [24], to solve compressible multiphase flows is also found in the literature [25,26,55,56], where numerical dissipation is largely reduced, and the moving interfaces,as well as other flow structures, can be resolved withsignificantly improvedaccuracy. However, implementinghigh orderschemes might generate numerical oscillations for compressible multiphase flows with complex EOS as discussed in [25], where even though the reconstructionswere carriedout in termsofthe characteristicvariables to reduce numericaloscillation,spu- rious disturbances are still observed when waves are reflected from the material interface. In a more recentwork [27], an intermediate state was introducedateach cell edgeincharacteristic decompositiontosuppressnumericaloscillations andstabilizecomputation.Furthermore,highordermonotonicity-preserving scheme[28] wasusedtoensurethebounded value for volume fraction. It is notedthat numerical dissipation even reduced in WENOand other high-order schemes stillremainsinconventionalEulerianadvectionschemes,whichmightbeproblematicinlong-termsimulations.Ingeneral, the implementationofhighorder shockcapturingschemesto interfacial compressiblemultiphaseflows demands further investigations.

There are different numerical methods to identify and compute moving interfaces incompressible multiphase flows, such as [29–31] on moving meshes and [32–36] on fixed meshes. As aforementioned, the VOF-type methods that use

1 Moreprecisely,itshouldbecalledsingle-statemodelorsingle-equivalent-fluid(SEF)model[21].WecallsuchmodelSEFinthepresentpaper.

thevolume fractionormassfractionastheidentificationfunction ofmovinginterface havebeenpopularlyused aswell, which are referred to as interface-capturing methods in our context. Interface-capturing methods resolve the interfaces onfixed Euleriangridsanduseadvectionschemestotransport thevolume/massfractionfunctions. Itiswell knownthat conventionalEulerianadvectionschemeshaveintrinsicnumericaldissipation andtendtosmearout thejumpsinvolume fractionormassfractionfunctionswhichare usedtoidentifythematerial interfacesbetweendifferentfluids. Inorderto keepthecompactthicknessofmaterialinterfacesduringcomputation,specialnumericaltechniquesareneededtosteepen thejumpsinthevolumeormassfractionfields.Forexample,in[37–39] theadvectionequationoftheinterfacefunctionis treatedbyartificialcompressionmethod.Asapost-processingapproach,anti-diffusiontechniqueshavealsobeenintroduced in [40] and [41]. An alternative approach is to reconstruct the volume fraction under the finite volume framework by usingspecialfunctions.TheTHINC(TangentofHyperbolaforINterfaceCapturing)method,forexample,usesthehyperbolic tangentfunction [42] tocapturethejumpsin volumefraction. Byvirtueofthedesirablecharacteristics ofthehyperbolic tangent function in mimicking the jump-like profile of the volume fraction field, the sharp interface can be accurately capturedinasimpleandefficientway[43,27]. However,when applyinginterface-sharping methodsto theSEFmodelsof multiphasecompressibleflows,velocityandpressureoscillationsmayoccuracrosstheinterface[39,43,41,38,44] duetothe inconsistencybetweenthephysicalvariablesandthevolumefractionfieldwithsharpenedorcompressedjumps.Asstated in[38,39],incontrasttoincompressibleflowswherethedensitiesoffluidsarefixed,artificialinterface-sharpeningscheme cannot be applied alone to volume fraction function in compressible multiphase flows. Modifications to other physical variables have to be made to maintain the consistency among the sharpened volume fraction and other physical fields [38,39,41,44].In [43],a homogeneousreconstruction hasbeenproposed wherethereconstructed volumefractionisused to extrapolate the remaining conservative variables across the interface to ensure the mechanical consistencyacross the isolatedmaterialinterfaces.

Thisworkpresentsanovelmethodologytoresolveproblem(II)addressedabove.Toalleviatethenumericaldissipation inhigh-resolution schemes forshockcapturing,a newschemeforspatial reconstruction hasbeen devisedto reduce nu- mericaldissipationsoastomaintainthesharpnessofthejumpsinvolumefractionthatidentifythemovinginterfaces.The scheme,so-calledMUSCL-THINC-BVD,implements theunderlying ideaofthe boundaryvariation diminishing (BVD)algo- rithm [45,46] with thetraditionalMUSCL(MonotoneUpstream-centered SchemesforConservation Laws)schemeandthe interface-sharpeningTHINCschemeastwobuilding-blocksforreconstruction.TheBVD algorithmchoosesareconstruction function betweenMUSCL andTHINC,so that the variations (jumps)of thereconstructed variables atcell boundariesare minimized,whicheffectivelyremovesthenumericaldissipationsinnumericalsolutions.Moresubstantially,MUSCL-THINC- BVDschemeisappliedtothevolumefractionandotherstatevariablessimultaneouslyinafinitevolumeframework,sothat theconsistencycanberealizedamongthephysicalvariablesandvolumefractionthroughoutthesolutionprocedure.Hence themanipulationsrequiredinotherexistingmethodstoadjustthephysicalvariablestobeconsistentwiththevolumefrac- tion,arenotneededinthepresentmethod.Thenumericalmodelisformulatedunderastandardfinitevolumeframework withaRiemannsolverinthewave propagationform[47].Thenumericaltestsverifythecapabilityofthepresentmethod incapturingthematerialinterfaceasawell-resolvedsharpjumpinvolumefraction.Thenumericalresultsforawiderange ofbenchmarktestsinone,twoandthreedimensionsshowsuperiorsolutionqualitycompetitivetootherexistingmethods.

Thispaperisorganizedasfollows.InSection 2,thegoverningequationsofthefive-equationmodelfortwo-phaseflow withmovinginterfacesaredescribed.InSection3,afterabriefreviewofthefinitevolumemethodwiththeRiemannsolver inwave-propagationformforsolvingthequasi-conservativefive-equationmodel,thedetailsofthenewMUSCL-THINC-BVD scheme forspatial reconstruction are presented.Numerical resultsof benchmarktestsare presented incomparison with otherhigh-resolutionmethodsinsection4.SomeconcludingremarksendthepaperinSection5.

2. Mathematicalmodel 2.1. Governingequations

In this work, the inviscid compressible two-component flows are formulated by the five-equation model developed in [21]. Byassuming that thematerialinterface isinequilibriumofmixedpressureandvelocity, thefive-equationmodel consistsoftwocontinuityequationsforphasicmass,amomentumequation,anenergyequationandanadvectionequation ofvolumefractionasfollows

∂

∂

t

( α

1

ρ

1

) + ∇ · ( α

1

ρ

1u

) =

0

,

∂

∂

t

( α

₂

ρ

₂

) + ∇ · ( α

₂

ρ

₂u

) =

⁰

,

∂

∂

t

( ρ

u

) + ∇ · ( ρ

u

⊗

^u

) + ∇

^p

=

⁰

,

∂

E

∂

t

+ ∇ · (

^Eu

+

^pu

) =

⁰

,

∂ α

1

∂

t

+

^u

· ∇ α

1

=

⁰

,

(1)

where

ρ

k and

α

k∈ [⁰,1] ^{denotein turnthekth phasic densityandvolumefractionfork}=¹,2,u the vector ofparticle velocity, p themixturepressureandE thetotalenergy.Whenconsideringmorethantwo-phases,thefive-equationmodel can beextendedby supplementingadditionalcontinuityequationsandvolumefractionadvectionequationsforeachnew phase.

2.2. Closurestrategy

Toclosethesystem,thefluidofeachphaseisassumedtosatisfytheMie–Grüneisenequationofstate, pk

( ρ

k

,

ek

) =

p_∞,k

( ρ

k

) + ρ

k



k

( ρ

k

) 

ek

−

e_∞,k

( ρ

k

) 

,

(2)

wherek= (¹/

ρ

k)(∂pk/∂ek)|ρk istheGrüneisencoefficient,and p_∞,k,e_∞,k aretheproperlychosen statesofthepressure and internal energyalong some referencecurves(e.g., alongan isentrope orother empirically fittingcurves) inorder to matchtheexperimentaldataoftheexaminedmaterial[48].Usually,parametersk,p_∞,kande_∞,kcanbetakenasfunctions only of the density. Thisequation of state can be employed to approximatea wide variety ofmaterials includingsome gaseousorsolidexplosivesandsolidmetalsunderhighpressure.

Theconservativenessconstraintsleadtothemixingformulaforvolumefraction,densityandinternalenergyasfollows,

α

1

+ α

2

=

¹

, α

1

ρ

1

+ α

2

ρ

2

= ρ , α

1

ρ

1e1

+ α

2

ρ

2e2

= ρ

e

.

(3)

Derivedin[20],themixtureGrüneisencoefficient,pressurep_∞ andinternalenergye_∞ canbeexpressedas

α

1



1

( ρ

1

) + α

2



2

( ρ

2

) =

¹

 ,

α

1

ρ

1e_∞,1

( ρ

1

) + α

2

ρ

2e_∞,2

( ρ

2

) = ρ

e_∞

, α

1

p_∞,1

( ρ

1

)



1

( ρ

1

) + α

2

p_∞,2

( ρ

2

)



2

( ρ

2

) =

^p∞

( ρ )

( ρ ) ,

(4)

undertheisobaricassumption.Themixturepressureisthencalculatedby

p

=

 ρ

e

−



2 k=¹

α

_k

ρ

_ke_∞,k

( ρ

_k

) +



2 k=¹

α

_k^p∞,k

⁽ ρ

_k

)



_k

( ρ

_k

)

  

2

k=¹

α

_k



_k

( ρ

_k

) .

(5)

ItshouldbenotedthatthemixingrulesofEq. (4) andEq. (5) ensurethatthemixedpressureisfreeofspuriousoscillations across the materialinterfaces [14,18,16,20,19]. Followingthe fiveequations modelunderisobaric closure[21], thesound speedofmixturecouldbecalculatedasthevolumetricaverageofthephasicsoundspeedsas

c²

= α

₁c²₁

+ α

₂c²₂

.

(6)

3. Numericalmethods

Forthesakeofsimplicity,thenumericalmethodinonedimensionisintroducedinthissection.Ournumericalmethod can be extended to multidimensions on structured grids directly in the dimension-wise reconstruction fashion. We will first review the finitevolume methodin thewave propagationform [47] usedin thiswork andthengive details of the newMUSCL-THINC-BVDreconstruction scheme.Comparedtothe classicalfinite-volumemethod,themodelsystemwhich is a quasi-conservativesystemof equationscan be approximatedin aconsistent andaccurate mannerwith wave propa- gationmethod[49,50,43].However,beingaspatialreconstructionschemeforhyperbolicfluxes,MUSCL-THINC-BVDcanbe implementedinanyfinitevolumeframeworkstraightforwardly.

3.1. Wavepropagationmethod

Werewritetheonedimensionalquasi-conservativefive-equationmodel(1) as

∂

q

∂

t

+ ∂

f

(

q

)

∂

x

+

^B

(

q

) ∂

q

∂

x

=

⁰

,

(7)

wherethevectorsofphysicalvariablesq andfluxfunctionsf are

q

= ( α

1

ρ

1

, α

2

ρ

2

, ρ

u

,

E

, α

1

)

^T

,

f

= ( α

1

ρ

1u

, α

2

ρ

2u

, ρ

uu

+

^p

,

Eu

+

^pu

,

0

)

^T

,

⁽⁸⁾

respectively.ThematrixB isdefinedas

B

=

diag

(

0

,

0

,

0

,

0

,

u

) ,

(9)

whereu denotesthevelocitycomponentinx direction.

WedividethecomputationaldomainintoN non-overlappingcellelements, Ci:^x∈ [^xi−¹/2,xi+¹/2]^,i=¹,2,. . . ,N,using auniformgridwiththespacingx=^xi+¹/2−^xi−¹/2.Forastandardfinitevolumemethod,thevolume-integratedaverage valueq¯_i(t)incellC_iisdefinedas

¯

q_i

(

t

) ≈

¹



x

x



_i+1/2

x_i−1/2

q

(

x

,

t

)

dx

.

(10)

Denoting all thespatial discretization termsin(7) by L( ¯^q(t)), thesemi-discrete version ofthe finitevolume formulation canbeexpressedasasystemofordinarydifferentialequations(ODEs)

∂

q

¯ (

t

)

∂

t

= L (¯

^q

(

t

)) .

(11)

Inthewave-propagationmethod,thespatialdiscretizationforcellCi iscomputedby

L ( ¯

qi

(

t

)) = −

¹



x

 A

⁺



qi−¹/2

+ A

⁻



qi+¹/2

+ A

^qi



(12)

whereA⁺q_i−1/2andA⁻q_i+1/2aretheright- andleft-movingfluctuations,respectively,andA^qiisthetotalfluctuation within Ci. Riemann problemsare solved to determinethese fluctuations. The right- and left-moving fluctuationscan be calculatedby

A

^±



q_i−1/2

=



3 k=1

s^k

q_i^L_−1/2

,

q_i^R_−1/2



_±

W

^k

q^L_i−1/2

,

q^R_i−1/2



,

(13)

wheremovingspeedss^kandthejumpsW^{k (k}=1,2,3)ofthreepropagatingdiscontinuitiescanbesolvedbytheRiemann solver[51] withthereconstructedvaluesq_i^L_−1/2 andq_i^R_−1/2 computedfromthereconstruction functions^q˜i−1(x)and^q˜i(x) totheleftandrightsidesofcelledgexi−¹/2,respectively.Similarly,thetotalfluctuationcanbedeterminedby

A

q_i

=



3 k=¹

s^k

q^R_i−1/2

,

q^L_i+1/2



_±

W

^k

q_i^R_−1/2

,

q^L_i+1/2



.

(14)

Wewilldescribewithdetailsaboutthereconstructionstogetthesevalues,q^L_i−1/2andq^R_i−1/2,atcellboundariesinthenext subsectionasthecorepartofthispaper.

In practice, given the reconstructed values q^L_i−1/2 and q^R_i−1/2, the minimum and maximum moving speeds s¹(q^L_i−1/2,q_i^R_−1/2)ands³(q_i^L_−1/2,q^R_i−1/2)arecomputedas

s¹

=

^min

{

^u_i^L_−1/2

−

^c_i^L_−1/2

,

u_i^R_−1/2

−

^c_i^R_−1/2

},

s³

=

^max

{

^u_i^L_−1/2

+

^c_i^L_−1/2

,

u_i^R_−1/2

+

^c_i^R_−1/2

},

⁽¹⁵⁾

wherec^L_i−1/2 andc^R_i−1/2 are thesoundspeedscalculatedby reconstructedvaluesq^L_i−1/2 andq_i^R_−1/2 respectively. Thenthe speedofthemiddlewaveisestimatedbytheHLLCRiemannsolver[51] as

s²

=

^p

R

i−¹/2

−

^pL_i−1/2

+ ρ

_i^L_−1/2u^L_i−1/2

(

s¹

−

^uL_i−1/2

) − ρ

_i^R_−1/2u^R_i−1/2

(

s³

−

^u_i^R_−1/2

)

ρ

_i^L_−1/2

(

s¹

−

u^L_i−1/2

) − ρ

_i^R_−1/2

(

s³

−

u_i^R_−1/2

) .

(16) Theleft-sideintermediatestatevariablesq^∗_i−^L_1/2isevaluatedby

q^∗_i−^L_1/2

= (

u_i^L_−1/2

−

^s1

)

q_i^L_−1/2

+ (

^pL_i−1/2^nL_i−1/2

−

^p∗_i−1/2^n∗_i−1/2

)

s²

−

^s1

,

(17)

where thevector u_i^L_−1/2= (^uLi−¹/2,u_i^L_−1/2,u^L_i−1/2,u^L_i−1/2,s²),n_i^L_−1/2= (⁰,0,1,u^L_i−1/2,0),n^∗_i−1/2= (⁰,0,1,s²,0)andthein- termediatepressuremaybeestimatedas

p^∗_i−1/2

= ρ

_i^L_−1/2

(

u^L_i−1/2

−

^s1

)(

u^L_i−1/2

−

^s2

) +

^pL_i−1/2

= ρ

_i^R_−1/2

(

u^R_i−1/2

−

^s1

)(

u_i^R_−1/2

−

^s2

) +

^p_i^R_−1/2

.

(18) Analogously,theright-sideintermediatestatevariablesq^∗_i−^R_1/2 is

q^∗_i−^R_1/2

= (

u^R_i−1/2

−

^s3

)

q^R_i−1/2

+ (

^p_i^R_−1/2ⁿ_i^R_−1/2

−

^p∗_i−1/2^n∗_i−1/2

)

s²

−

^s3

.

(19)

ThenwecalculatethejumpsW^k(q^L_i−1/2,q_i^R_−1/2)as

W

¹

=

^q∗i−^L1/2

−

^qLi−¹/2

, W

²

=

q^∗_i−^R_1/2

−

q^∗_i−^L_1/2

, W

³

=

q^R_i−1/2

−

q^∗_i−^R_1/2

.

(20)

Forthenon-conservativeterm,thejumpforthevolumefractionissimplyzeroforbothW^1,W^3and

α

₁^R−

α

₁^L forW^2from above equations.Giventhe spatial discretization, we employ thethree-stage third-order SSP (Strong Stability-Preserving) Runge–Kuttascheme[52]

¯

q^∗

= ¯

^qn

+ 

^t

L 

¯

qⁿ



,

¯

q^∗∗

=

³

4q

¯

ⁿ

+

¹ 4q

¯

^∗

+

¹

4



t

L 

¯

q^∗



,

¯

qⁿ⁺¹

=

¹

3q

¯

ⁿ

+

² 3q

¯

^∗

+

²

3



t

L 

¯

q^∗∗



,

(21)

to solve thetime evolution ODEs,where q¯^∗ andq¯^∗∗ denotethe intermediate valuesatthe sub-steps. Fortemporalinte- gration, we havealso tested the Eulerfirst-order explicit andsecond order Runge–Kutta schemes,which produce stable numerical results as the third order SSP Runge–Kutta does. It is observed that the results from first andsecond-order temporalschemesareslightlydiffusive.

3.2. MUSCL-THINC-BVDreconstruction

Intheprevioussubsection,theboundaryvalues,q_i^L_−1/2 andq_i^R_−1/2,arelefttobedetermined,whichwillbepresentedin thissubsection.Wedenoteanysinglevariableforreconstructionbyq,whichcanbeprimitivevariable,conservativevariable orcharacteristicvariable.

Thevaluesq^L_i−1/2andq_i^R_+1/2atcellboundariesarecomputedfromthepiecewisereconstructionfunctionsq˜i(x)incell Ci. In the presentwork, theMUSCL-THINC-BVD reconstruction schemeis designedto captureboth smooth andnon-smooth solutions.TheBVDalgorithmmakesuseoftheMUSCLscheme[23] andtheTHINCscheme[42] asthecandidatesforspatial reconstruction.

In theMUSCLscheme,a piecewiselinear functionisconstructed fromthe volume-integratedaveragevaluesq¯i,which reads

˜

q_i

(

x

)

^{MU SC L}

= ¯

^qi

+ σ

i

(

x

−

^xi

),

(22)

where x∈ [^xi−1/2,x_i+1/2] ^and

σ

i isthe slope definedatthe cellcenter x_i= ¹₂(x_i−1/2+^xi+1/2).Toprevent numerical os- cillation, aslope limiter [23,50] isused to getnumerical solutionssatisfying theTVD property.The reconstructed values atcellboundariesfromMUSCLreconstructionaredenotedasq_i^R₋^{,MU SC L}_1/2 = ˜^qi(xi−¹/2)^{MU SC L} andq^L_i+^{,MU SC L}_1/2 = ˜^qi(xi+¹/2)^{MU SC L}. The MUSCLscheme,inspiteofpopularuseinvariousnumericalmodels,hasexcessivenumericaldissipationandtendsto smearoutflowstructures,whichmightbeafataldrawbackinsimulatinginterfacialmultiphaseflows.

Being another reconstruction candidate, the THINC scheme [42,53] uses the hyperbolic tangent function, which is a differentiable andmonotonefunction thatfitswell astep-like discontinuity.TheTHINC reconstructionfunctionis written as

˜

q_i

(

x

)

^{T H I NC}

= ¯

^qmin

+

q

¯

_max 2

1

+ θ

^tanh

β

x

−

xi−¹/2

x_i+1/2

−

^xi−1/2

− ˜

^xi



,

(23)

where q¯min=^min(¯qi−¹,q¯i+¹),q¯max=^max(¯qi−¹,q¯i+¹)− ¯^qmin andθ=^sgn(¯qi+¹− ¯^qi−¹).The jumpthicknessiscontrolledby parameter β.Same as[45], β=¹.6 performs well for all numericaltestspresented inthis paper.When β becomes too smaller, theresultstendto bemorediffusive whilealarger β enforcesthe anti-diffusioneffectandtendsto steepenthe jumps in the numerical solutions. From our numerical experiments, a β valuedfrom 1.4 to 2.0 isable to give good or acceptable results. In our numerical tests shown later a constant value of β is used in one dimensional problems. For

Fig. 1. Illustration of one possible situation corresponding to|^qLi−^{,MU SC L}1/2 −^qi^R−^{,T H I NC}1/2 | + |^qi^L+^{,T H I NC}1/2 −^qi^R+^{,MU SC L}1/2 |when calculating T B V_i^{T H I NC}_,min .

multidimensionalcases,β dependsonthediscontinuityorientation.Followingtheworkin[53],parameterinx directionis computedbyβx= β|ⁿx|^wherenx isthex componentoftheunit normalofthediscontinuity.Likewise,theparametersin y andz directions,βyandβz,canbedeterminedwiththecorrespondingcomponentsoftheunitnormalvectorny andnz inasimilarway.Inordertokeepsimplicity,atraditionalYoung’s algorithmdescribedin[54] isusedtocalculatetheunit normal.Theunknown ˜^xi,whichrepresentsthelocationofthejumpcenter,iscomputedfrom^q¯i=_¹_x_xi+1/2

xi−1/2 ^q˜i(x)^{T H I NC}dx.

ThenthereconstructedvaluesatcellboundariesbyTHINCfunctioncanbeexpressedby q^L_i+^{,T H I NC}_1/2

= ˜

qi

(

xi+¹/2

)

^{T H I NC}

= ¯

qmin

+

q

¯

max

2

1

+ θ

^tanh

(β) +

A 1

+

^{A tanh}

(β)

 ,

q^R_i−^{,T H I NC}_1/2

= ˜

qi

(

xi−¹/2

)

^{T H I NC}

= ¯

qmin

+

q

¯

max

2

(

1

+ θ

A

) ,

(24)

where A= ^B/cosh_tanh^(β)−_(β) ^1,B=^exp(θ β(2 C−¹)),andC=^q¯i− ¯^qmin+

¯

qmax+

^{isamappingfactortoprojectthephysicalfieldsonto} [⁰,1]^{.Asmallpositive,} =^{10−20isintroducedtopreventarithmeticfailure.}

The final reconstruction function isdetermined by the BVD algorithm, whichchooses the reconstruction function be- tweenq˜i(x)^{MU SC L} andq˜i(x)^{T H I NC} sothatthevariations ofthereconstructedvaluesatcellboundariesareminimized. BVD algorithm prefers the THINC reconstruction ^q˜i(x)^{T H I NC} within a cell wherea discontinuityexists. It is sensiblethat the THINCreconstructionshouldonlybeemployedwhenadiscontinuityisdetected.Inpractice,wemakeuseofthefollowing conditionsasanadditionalcriteriontoimplementtheTHINCreconstruction

δ <

C

<

1

− δ

^and

(¯

q_i+1

− ¯

^qi

)(

q

¯

_i

− ¯

^qi−1

) >

0

,

(25)

whereδisasmallpositive.

Inpresentwork,weuseamodifiedBVDalgorithmthatdeterminesthereconstructionfunctionby

˜

qi

(

x

)

^{B V D}

=



^q

˜

i

(

x

)

^{T H I NC} if

δ <

C

<

1

− δ,

^and

(

^q

¯

i+¹

− ¯

^qi

)(¯

^qi

− ¯

^qi−¹

) >

0

,

and T B V_i^{T H I NC}_,min

<

T B V_i^{MU SC L}_,min

˜

q_i

(

x

)

^{MU SC L} otherwise

,

(26)

wheretheminimumvalueoftotalboundaryvariation(TBV)T B V_i^P_,minforreconstructionfunctionq˜i(x)^P iscomputedby

T B V_i^P_,min

=

^min

(|

^qL_i−^{,MU SC L}_1/2

−

^q_i^R₋^,P_1/2

| + |

^qL_i+^,P_1/2

−

^qR_i+^{,MU SC L}_1/2

|, |

^qL_i−^{,T H I NC}_1/2

−

^q_i^R₋^,P_1/2

| + |

^q_i^L₊^,P_1/2

−

^q_i^R₊^{,T H I NC}_1/2

|,

|

q^L_i−^{,MU SC L}_1/2

−

q_i^R₋^,P_1/2

| + |

q^L_i+^,P_1/2

−

q^R_i+^{,T H I NC}_1/2

|, |

q^L_i−^{,T H I NC}_1/2

−

q^R_i−^,P_1/2

| + |

q^L_i+^,P_1/2

−

q_i^R₊^{,MU SC L}_1/2

|),

⁽²⁷⁾

where P standsforeitherT H I NC orMU SC L.

Hence,THINCreconstructionfunctionwillbeemployedinthetargetedcelliftheminimumTBVvalueofTHINCissmaller than that of MUSCL. In Fig.1, we illustrate one possible situation corresponding to |^qL_i−^{,MU SC L}_1/2 −^qR_i−^{,T H I NC}_1/2 |+ |^qL_i+^{,T H I NC}_1/2 − q^R_i+^{,MU SC L}_1/2 | ^{whenevaluating T B V}i^{T H I NC},min . Asstatedin [45],the BVD algorithmwill realizethepolynomial interpolation for smooth solution whilefordiscontinuous solution astep like functionwill be preferred. Itis notedthat thepresentBVD algorithm, (26) and (27), is slightly different from that in [45]. The present BVD algorithm minimizes the total BVs at two ends of the target cell. Ournumerical tests show that the presentBVD algorithm can effectively reduce numerical dissipationsandpreventtheflowstructuresfrombeingsmearedoutasthatin [45].

Fig. 2. NumericalerrorsindensitycomputedbythepureMUSCLscheme(a)andtheMUSCL-THINC-BVDscheme(b).Solidanddashedlinesindicatethe slopesof−^{1 and}−2 respectively.

We have evaluated the numerical errors and convergence rate of the MUSCL-THINC-BVD scheme using the mesh- refinementtests.Propagationofacousticwavesinasinglegasisconsidered.Theinitialconditionissetasthesameas [26], whereaninitialperturbationisaddedtodensityandpressurefieldby

ρ

₀

₌

1

+ η

h

(

x

),

u₀

=

⁰

,

p₀

=

¹

γ ⁺ η

h

(

x

),

(28)

wheretheperturbationfunction isspecifiedash(x)=^sin8(

π

x)and

η

=^10−4.Theconvergencerateisstudiedby refining the meshfromN=^{10 to N}=^{80 andevaluated forthedensityfield.The L}1 errorsforbothMUSCL schemeandMUSCL- THINC-BVD schemeare summarized inFig. 2.As expected, the MUSCL-THINC-BVD schemerealizes the convergencerate of thepolynomial-based reconstructioncandidate forsmoothsolutions, andshowsaconvergence behaviorsimilar tothe MUSCL scheme.Asdiscussed anddemonstratedin[45],theBVD algorithmprefersthehigh-order polynomialreconstruc- tion forsmooth solution,andusing a high-orderpolynomial reconstruction, such asthe WENOscheme in[45],leads to high-order convergencerate. Eventhough theMUSCL schemeisused inthepresentBVD algorithm andtheconvergence rate isthus under 2nd order, thenumerical errors inthe resultsof MUSCL-THINC-BVD schemeare significantly reduced comparedtothepureMUSCLscheme.Shownin[45],aswellasthenumericalresultsformultiphaseflowsinthispaper,by minimizingthereconstructeddifferencesacrosscellboundaries,aBVDschemecaneffectivelyreducenumericaldissipation andthusbetterresolvevorticalstructures,suchastheKelvin–Helmholtzassociatedwithdensitydiscontinuities.

It isnotedthat the MUSCL-THINC-BVDschemecan be usedforsingle-phaseflows asshowcasedin [45] andgetsub- stantiallyimprovedsolutionquality. FortheSEFmodelofmultiphaseflows,ifamixing/averagingmodel,suchasthosein [14–21], isdesignedtomaintain thephysicalbalanceacrossmaterialinterfaces,implementingMUSCL-THINC-BVDscheme totheresultingSEFmodelwillnotgeneratespuriousoscillationsinthenumericalsolution.

Asshownlaterinthenumericaltests,discontinuitiesincludingmaterialinterfacecanberesolvedbytheMUSCL-THINC- BVD scheme with substantially reduced numerical dissipation in comparison with other existing methods. The material interfacecanbecapturedsharply,whileanyextrastep,likeanti-diffusionorotherartificialinterfacesharpeningtechniques used intheexistingworks[27,43,38,39], isnotneededhere. Moreimportantly,theMUSCL-THINC-BVDschemeisapplied not onlytothevolumefractionbutalsoto otherphysicalvariables,suchasthephasicdensity,whichautomaticallyleads totheconsistencyamongthereconstructedphysicalfields.Asobservedinournumericalresults,nospuriousnumericalos- cillationisgeneratedinthevicinityofmaterialinterfaces.Itisusuallynottrivialtopreventnumericaloscillationsforother anti-diffusion or artificialcompressionmethods aforementioned.For example,in[44,39,43] anti-diffusion post-processing stepsarerequiredtoadjustthestatevariablesacrossthematerialinterfacestogetaroundtheoscillations.

4. Numericalresults

Numericaltestsinone-,two- andthree-dimensionsarepresentedinthissectiontoverifytheproposedMUSCL-THINC- BVDschemeincomparisonwiththeWENOscheme.HereweusetheWENOschemein[24] whichisoneofrepresentative highordershock-capturingschemes.WedenoteitasWENO-JSinourtests.Asaddressedin[25],theWENOreconstruction shouldbe implementedforcharacteristicfieldsinordertoreducethenumericaloscillations,whichisnot aneasytaskfor