A conservative interface method for compressible flows

www.elsevier.com/locate/jcp

Aconservativeinterfacemethodforcompressibleflows

X.Y.Hu

abca,*,B.C.Khoob,N.A.Adamsa,F.L.Huang

cLehrstuhlfu¨rAerodynamik,TechnischeUniversita¨tMu¨nchen,85747Garching,Germany

DepartmentofMechanicalEngineering,NationalUniversityofSingapore,117576Singapore,Singapore

NationalKeyLaboratoryforPreventionandControlofExplosiveDisasters,BeijingInstituteofTechnology,100080Beijing,China

Received21October2005;receivedinrevisedform22February2006;accepted3April2006

Availableonline19May2006

Abstract

Inthiswork,wepresentaconservativeinterfacemethodforbothmulti-fluidandcomplexboundaryproblems,inwhichthestandardfinitevolumeschemeonCartesiangridsismodifiedbyconsideringcomputationalcellsbeingcutbyinterface.Whilethediscretizedgoverningequationsareupdatedconservatively,themethodtreatsthetopologicalchangesnaturallybycombininginterfacedescriptionandgeometricoperationswithalevelsettechnique.Extensivetestsin1Darecarriedout,and2Dexamplessuggestthatthepresentschemeisabletohandlemulti-fluidandcomplex(staticormoving)bound-aryproblemsinastraightforwardwaywithgoodrobustnessandaccuracy.Ó2006ElsevierInc.Allrightsreserved.

Keywords:Compressibleflow;Multi-fluidproblem;Complexboundary;Levelsetfunction

1.Introduction

Variousnumericalmethodshavebeendevelopedtosimulateandstudythedynamicsofcompressibleflowsathighspeeds.Inallcases,asignificantdifficultyencounteredwiththesenumericalcomputationsisthetreat-mentofamaterialinterface.Ingeneral,therearetwomaintypesofinterfacerelatedproblems:oneisthemulti-fluidprobleminwhichtheinterfaceseparatesinteractingdifferentfluids;theotheristhecomplexboundaryprobleminwhichtheinterfacedefinesacomplex,staticormoving,boundarybetweenthefluidandarigidsolidwallorabody.

Usually,theinterfacesforthesetwotypesofproblemsaretreatedinquitedifferentways.Foramulti-fluidproblem,theinterfaceisoftendefinedonaCartesiangridasatransitionregionwithasteepgradient.Somenumericalmethodsusingthis‘smeared’interfacerepresentationarethevolumeoffluid(VOF)method[1,18,49]orconstrainedinterpolationprofile(CIP)method[46].Foracomplexboundaryproblem,theinter-faceisoftendefinedasanon-smearedboundaryofastructuredbody-fittedgrid[4],oranunstructuredgrid[23,28,36].Ontheotherhand,therearestillalternativeapproachesinwhichtheabovetwotypesofinterfaces

*Correspondingauthor.Tel.:+49216152;fax:+49216139.E-mailaddress:xiangyu.hu@aer.mw.tum.de(X.Y.Hu).

0021-9991/$-seefrontmatterÓ2006ElsevierInc.Allrightsreserved.doi:10.1016/j.jcp.2006.04.001

554X.Y.Huetal./JournalofComputationalPhysics219(2006)553–578

arebothtreatedonaCartesiangrid.Oneapproachistoreconstructa‘‘smeared’’complexboundarywiththeVOFtechnique[19].Theotherapproachistotrackbothmulti-fluidinterfacesandcomplexboundarieswithanon-smearingrepresentation.Glimmetal.[15]trackthemulti-fluidinterfaceasanon-smearedinternalmov-ingboundarywithagridintersectiontechnique.Similargeometrictreatmentsoftheboundaryarealsodevel-opedforCartesianapproachestocomplexboundaryproblems[5,7,38,47,48].Eventhoughthesementionedalternativeapproachesareconceptuallysimpletherearetwomajorlimitations:oneisthatdefiningandtrack-ingoftheinterfacerequirerathercomplicatedprocedureswiththegridintersectiontechnique,especiallyinthecasesoflargetopologyvariations;theotheristhelackofconservativeproperties,whichusuallyleadstolowaccuracyneartheinterface,especiallywhenstrongorlongtimescaleinterfaceinteractionsareinvolved.Thelevelsettechnique[10,34,43]isareasonablyeasywaytodefineandtrackanon-smearedmulti-fluidinterface.Todealwiththeissueofconservation,theoriginalworksofBergerandLeVque[5]andQuirk[38]oncomplexboundaryproblemsdonottreattheinterfacespecificallybutuseacomplexadaptivemeshrefinementstrategytoincreaseaccuracyneartheboundary.Otherapproaches,e.g.GFM(ghostfluidmethod)andrelatedmethods[8,12,20,27]employspecialinternalboundaryconditionsatthemulti-fluidinterfaceintheattempttocorrectorreducedconservationerrors.Fairlysimilarprocedurescanalsobeimplementedforcom-plexboundaryproblems[2,13,14,21].Whiletheseapproacheshaveledtopromisingresultstheypossessstrictlyonlyafirst-orderconvergenceratefortheconservationerrors.TohandlethestabilitydifficultyofaconservativediscretizationwiththeVOFmethod,Colella[9]andMillerandColella[32]introducedaredis-tributionapproachformulti-fluidproblemsinwhichconservationdefectsresultingfromnon-conservativeupdatesareredistributedontoneighboringcellsaccordingtomassweighting.Falcovitzetal.[11]satisfycon-servationnearacomplexmovingboundarybyaStrang-typeoperatorsplitting.Glimmetal.[16]introducedaconservativemethodbasedonagridintersectiontechnique,inwhicheachsmall-sizecutcellismergedwithoneofitsneighboringcells,formulti-fluidflowswiththeassumptionofnotopologicalchangeinthesolution.Inthispaper,weproposeaCartesianinterfacemethodsuitableforbothmulti-fluidandcomplexboundaryproblems.Weuseastandardfinitevolumeschemeforthefarinterfaceregionwhichisslightlymodifiedforthenearinterfaceregion.UnlikeMillerandColella[32],thepresentmethodupdatesthediscretizedgoverningequationsfullyconservativeforbothfluidsindividuallyandforinterfaceexchangesinmulti-fluidproblems.Smallcutcellsaretreatedwithamixingprocedurewhich,however,isdifferentfromtheapproachesfoundin[32,16].Asourmixingprocedureisseparatedfromthefluxupdatingofconservativevariables,itiseasytoimplement,especiallywhencomplexgeometriesareconsidered.Also,unliketheworkofFalcovitzetal.[11],thepresentmethodoffersafairlysimpleimplementationinmulti-dimensionandmulti-leveltimeintegra-tionswithoutspace–timesplitting.Theassumptionofnon-topologicalchangeasin[16]isnotnecessaryinourformulation.Furthermore,asthelevelsettechniqueisusedfortheinterfacedescriptionandgeometricoper-ations,themethodmaintainsthesimplicityofGFM-likemethods.Indeed,inourimplementation,weusetheinterfaceinteractionmethod[20]toobtaintheinterfaceconditionsforthenearinterfacegridpoints.However,theseinterfaceinteractionsarenotemployedfordefiningtheghostnodestatesbutforthedirectcalculationofmomentumandenergyexchangesacrosstheinterface.Accordingly,theghostnodeisfilledonlyforthestencilinterpolationratherthanforimplementingtheinternalboundaryconditioncommonforGFM-likemethods.2.Overviewofthemethod

Assumingthefluidisinviscidandcompressiblethegoverningequationoftheflowcanbewrittenasasys-temofconservationlaws

oU

þrÁF¼0ot

onX;

ð1Þ

whereUisthedensityoftheconservedquantitiesofmass,momentumandtotalenergy,andFrepresentsthecorrespondingfluxfunctions.WhenaninterfaceC(t)separatesthedomainXintotwosub-domainsX1(t)andX2(t),asforamulti-fluidproblem,theevolutionoftheinterfaceisdeterminedbytheinterfaceconditiongivenbyatwo-materialRiemannproblem

RðUfluid1;Ufluid2Þ¼0onCðtÞ:

ð2Þ

X.Y.Huetal./JournalofComputationalPhysics219(2006)553–578555

Whenacomplexboundaryproblemisconsidered,theevolutionoftheinterfaceisdeterminedbythebound-aryvelocityvrg.Thepressureontheinterfaceisthendeterminedbytheone-sidedRiemannproblem

RðUfluid;vrgÞ¼0

onCðtÞ:

ð3Þ

Theboundaryvelocitycanbeprescribed.Iftheassociatedrigidsolidbodymovesaccordingtoaninertialcou-plingtheaccelerationarg(t)canbecalculatedby

Z1

pdsð4ÞargðtÞ¼

mrgCðtÞinwhichtheintegrationiscarriedoutovertheentireinterface,andmrgisthesolidbodymass.

FollowingMillerandColella[32],weconsiderEq.(1)forthefluidoccupyingthesub-domainX1onatwo-dimensionalCartesiangridwithgridspacingsDxandDy.AfinitevolumediscretizationcanbeobtainedfromintegratingEq.(1)overthespace–timevolumeDij˙X1(t)ofacomputationalcell(i,j)occupiedbythefluid

󰀂󰀃Znþ1Z

oU

dtdxdyþrÁF¼0;ð5Þ

otnDij\\X1ðtÞwhereDij=DxDyisthecellvolume.Dij\\X1(t)canberepresentedbyai,j(t)DxDywhereai,j(t)isthetime

dependentvolumefractionoftheconsideredfluidandsatisfying1PaP0.ByanapplicationofGauss’sthe-orem,oneobtains

Znþ1ZZnþ1Z

oU

dtdxdydtdxdyFÁn¼0;ð6Þþ

otnnDij\\X1ðtÞoDij\\CðtÞwhereoDijarethefourcellfacesintersectingorthogonallywiththegridatfourlocations(xi+Dx/2,yj),

(xi,yj+Dy/2),(xiÀDx/2,yj)and(xi,yjÀDy/2).DenotingtheinterfacelocationasC(t),asshowninFig.1,oDij\\C(t)canberepresentedbytwoparts:oneisthecombinationofthefoursegmentsofthecellfacesbeingcutbytheinterface,whichcanbewrittenintheformofAi+1/2,j(t)Dy,Ai,j+1/2(t)Dx,AiÀ1/2,j(t)DyandAi,jÀ1/2(t)Dxwhere1PAP0istheaperture;theotheroneisthesegmentoftheinterfaceDCi,j(t)insideofthecell(i,j).Hence,Eq.(6)canberewrittenas

In(i,j+1)In+1Ai,j+1/2Ai-1/2,jAi+1/2,j=0(i,j)i,j(i-1,j)(i+1,j)Ai,j-1/2(i,j-1)ni,j+n+1i,jFig.1.Schematicofconservativediscretizationforacutcell.556X.Y.Huetal./JournalofComputationalPhysics219(2006)553–578

󰀄

þ1nþ1ani;jUi;jn

Àani;jUi;j

󰀅

i1h^^dtAiþ1=2;jðtÞFiþ1=2;jÀAiÀ1=2;jðtÞFiÀ1=2;jDt¼

Dxn

Znþ1iZnþ1Ã1h1^Â^i;jþ1=2ÀAi;jÀ1=2ðtÞF^i;jÀ1=2þþAi;jþ1=2ðtÞFXCi;jðtÞ;dtdt

DyDxDynn

nþ1

Z

ð7Þ

whereDtisthetimestepsize.ai,j(t)Ui,jandUi,jaretheconservativequantitiesinthecutcellandthecell-aver-^istheaveragecell-facefluxandageddensityofconservativequantitiesoftheconsideredfluid,respectively.F

^½Ci;jðtÞ󰀄istheaveragemomentumandenergyexchangeacrosstheinterfacesegmentdeterminedbytheinter-X

faceinteractionofEq.(2)or(3).Withexplicitfirst-orderforwardtimedifference,theaboveequationcanbeapproximatedas

iDt^Dthnnþ1nþ1n^iÀ1=2;jÀAiþ1=2;jF^iþ1=2;jXðDCi;jÞþAiÀ1=2;jFai;jUi;j¼ai;jUi;jþ

DxDyDxiDth^^þAi;jÀ1=2Fi;jÀ1=2ÀAi;jþ1=2Fi;jþ1=2:ð8Þ

DyNotethatallthetermsontheright-handsideareevaluatedattimestepn.Forfullcellswhicharenotcutbyaninterface,volumefractionsandaperturesbecomeunityandthecorrespondinginterfacesegmentsDCi,jvanish.Eq.(8)thensimplifiesto

󰀅Dt󰀄󰀅Dt󰀄^nþ1n^^^FiÀ1=2;jÀFiþ1=2;jþFi;jÀ1=2ÀFi;jþ1=2ð9ÞUi;j¼Ui;jþ

DxDyrecoveringastandardfinitevolumeschemeonatwo-dimensionalCartesiangrid.Ontheotherhand,since

beingeffectiveonlyincuttingcells,Eq.(8)canalsobeviewedasaslightmodificationofEq.(9)neartheinterface.

ItcanbefoundthatglobalconservationfortheconsideredfluidissatisfiedbysummingEq.(8)overtheentiredomain

X

i;j

þ1nþ1ani;jUi;j

¼

X

i;j

n

ani;jUi;j

þ

X

i;j

Dt^

XðDCi;jÞþboundaryterms;DxDy

ð10Þ

wheretheinterface-exchangeterm(secondtermonright-handside)vanishesinallfullcells.Notethatfortwo-fluidproblems,theevolutionofeachfluidiscomputedbysolvingEq.(8)wheretheinterface-exchangetermhasoppositesign.Therefore,overallconservationcanbeachievedbysummingEq.(10)forthetwofluids

XXXX

nnþ1nþ1

ai;jUi;j¼anð11Þi;jUi;jþboundaryterms:

k¼1;2

i;j

k¼1;2

i;j

Here,wemakethefollowingremarks:

󰀅AsstandardCartesianfinitevolumeschemesareusedforthediscretizationofdifferentialoperators,the

locationofthecentroidofacutcellisalwaysapproximatedbythelocationofagridpoint.Thisapproachhasbeenusedwithmuchsuccesspreviously[35,32].

󰀅Wheneverthevolumefractionofacutcellbecomessmallamixingprocedureisappliedtoavoidanextre-melysmalltimestep.

󰀅ForhigherordertimeintegrationaRunge–Kuttaschemecanbeemployedwhilemaintainingdiscretecon-servationsinceeveryRunge–Kuttasub-stepcanbeformulatedasEq.(8).

󰀅Theextensionoftheaboveformulationstothreedimensionsisfairlystraightforward.

Intheremainderofthepaperweelaborateontheaboveideasandimplementations.ThedescriptionofinterfacerelatedissuescanbefoundinSection3.Ourmethodisbasedonthelevelsettechnique.Besidestheinterfaceadvection,re-initializationandstateextrapolationfollowpreviousworks[12,20].Newmethods

X.Y.Huetal./JournalofComputationalPhysics219(2006)553–578557

areproposedtocalculatethecell-faceaperturesandvolumefractions.InSection4anovelbutsimpleandlowdissipationmixingprocedureisintroduced.InSection5theinterfaceinteractiontermsarederived.InSection6thedetailedimplementationstepsaregiveninthecontextofapplicationstomulti-fluidandcomplexbound-aryproblems,andpossibleextensiontothreedimensions.Finally,numericalexamplesinoneandtwodimen-sionsandbriefconcludingremarksaregiveninSections7and8,respectively.3.Interfacedescription3.1.Levelsetfunction

WeassociatethecomputationaldomainXwithasigneddistancefunction/(x,y,t),thatis$|/|=1,calledthelevelsetfunction[34].Knowing/wemaylocatetheinterfacebyfindingthezerolevelsetof/.ThatisC(t)={x,y:/(x,y,t)=0}whichdividestheentiredomainintotwosub-domains,eachofwhichcorrespondstoafluidorasolidbody,withoppositesignsofthelevelsetfunction/.Inthispaper,theapertureandvolumefractionforacomputationalcellaredefinedasA+anda+forthefluidregioncorrespondingto/>0,whileAÀandaÀforthefluidregioncorrespondingto/<0.

Forsimplegeometriesdescribedwithlinesorcircles,theinitiallevelsetfunctioncanbeconstructeddirectlybyusingoneormoredistancefunctionstothegivencurves.Ifmorecomplexinterfacesareconsideredtheinitiallevelsetfunctionisdifficulttobespecifiedanalytically,andusuallyneedstobeconstructednumericallyfromasetofconnectededgesorsurfacepatchesonvertices.Thefirststepfortheconstructionistodeterminewhetheragridpointisinsideoroutsideoftheinterface.Manypointclassifications,suchastherayintersec-tionmethod[22],canbeeasilyimplemented.Next,wefindapointontheinterfacehavingaminimumdistancetothegridpoint,anddesignatetheminimumdistanceandthedirectiontothispointasthemagnitudeof/anditsnormaldirection,respectively.

Withthelevelsetinitializedatthegridpointsitcanbefoundthatthecontinuousupdatingof/isequiv-alenttotheadvectionoftheinterfacebytheequation

/tþu/xþv/y¼0;

ð12Þ

whereuandvarethevelocitycomponentsforthelevelsetsinthexandydirections,respectively.High-orderschemes[12,40]havebeendevelopedtosolveEq.(12).Inpractice,thelevelsetsareonlyupdatedinthenearinterfaceregion,whichusuallyincludesthefirstandsecondnearestcell-layers.Fortheregionfarfromtheinterfacethelevelsetsarere-initialized[43]bytheequation

/sþsgnð/Þðjr/jÀ1Þ¼0;

ð13Þ

wheresgn(/)isasignfunction,tomaintainthesigneddistancepropertyoflevelsetfunction.Toallowforstencilinterpolationneartheinterfaceintermsofstandardfinitevolumeschemes[12]andsolvingfortheinterfaceconditionsongridpointsinnearinterfaceregion[20],thefluidstatesareextrapolatedtotheothersideoftheinterfacebytheextendingequation

qsÆNÁrq¼0:

ð14Þ

HereqistheextendedfluidstateandN”(Nx,Ny)isthenormaldirectionofthelevelset.Onmaynotethat+Nisusedtoextendquantitiesfromsub-domainwith/<0tosub-domainwith/>0,whileÀNisusedtoextendquantitiesfromsub-domainwith/>0tosub-domainwith/<0.Foragivenqtheextendingequationcanbesolveduntilasteadysolutioninthenearinterfaceregionisreached.3.2.Cell-faceapertures

Inthispaper,thecell-faceaperturesarecalculatedaccordingtothelevelsetdistributionalongthecellface.Onatwo-dimensionalCartesiangrid,thelevelsetfunctionatthecornersofacellcenteredat(i,j)isapprox-imatedby

558X.Y.Huetal./JournalofComputationalPhysics219(2006)553–578

Á1À

/i;jþ/iþ1;jþ/i;jþ1þ/iþ1;jþ1;4

Á1À

/iþ1=2;jÀ1=2¼/i;jþ/iþ1;jþ/i;jÀ1þ/iþ1;jÀ1;

4ð15Þ

Á1À

/iÀ1=2;jþ1=2¼/i;jþ/iÀ1;jþ/i;jþ1þ/iÀ1;jþ1;

4

Á1À

/iÀ1=2;jÀ1=2¼/i;jþ/iÀ1;jþ/i;jÀ1þ/iÀ1;jÀ1:

4

Asignchangeofthelevel-setfunctionalongacellfaceimpliesthatthecellfaceiscutbytheinterface.Forsuchcells,byassumingalineardistributionof/thecell-faceaperturescanbecalculatedas

/iþ1=2;jþ1=2¼Aiþ1=2;j¼aiþ1=2;jAiÀ1=2;j¼aiÀ1=2;jAi;jþ1=2¼bi;jþ1=2Ai;jÀ1=2¼bi;jÀ1=2

j/

if/iþ1=2;jþ1=2>0if/iÀ1=2;jþ1=2>0if/iþ1=2;jþ1=2>0if/iþ1=2;jÀ1=2>0

j

else1Àaiþ1=2;j;else1ÀaiÀ1=2;j;else1Àbi;jþ1=2;else1Àbi;jÀ1=2;

j/

j

ð16Þ

1=2;jþ1=21=2;jÆ1=2

whereaiÆ1=2;j¼j/iÆ1=2;jþ1iÆandbi;jÆ1=2¼j/iþ1=2;jÆ1iþ.Fig.2showstheschematicofthecalcula-=2jþj/iÆ1=2;jÀ1=2j=2jþj/iÀ1=2;jÆ1=2jtionforAi+1/2,j.Iflevelsetdoesnotchangesignalongacellface,theapertureiseither1(positivesign)or0(negativesign).Notethattheabovecell-faceapertures,denotedasA+,areforthefluidregioncorrespondingto/>0.Forthefluidcorrespondingto/<0,thecell-faceaperturesareAÀ=1ÀA+.3.3.Volumefractions

Withthestandardlevelsettechnique[43],thevolumefractionaþi;jforthefluidcorrespondingto/>0can

þ

beapproximatedbyasmoothedHeavisidefunctionai;j¼Hð/i;j;eÞwhereeisasmallpositivenumbersuchasthegridsize.ThevolumefractionaÀforthefluidcorrespondsto/<0canbecalculatedbytherelationaÀ=1Àa+.AsimpleformofthesmoothedHeavisidefunctioncanbewrittenas

(i,j+1)(i+1/2,j+1/2)Ai+1/2,j(i-1,j)(i,j)(i+1,j)(i+1/2,j-1/2)(i,j-1)Fig.2.ThecalculationoftheapertureAi+1/2,j.X.Y.Huetal./JournalofComputationalPhysics219(2006)553–578559

Hð/;eÞ¼

8><0

1>2:

/þ2þ21psineÀp/Á

e/<Àe;Àe6/6e;/>e;

ð17Þ

1

whichdependsonthelevelsetvaluesonly.Byincludingmoreinformationonthelevelsetnormaldirectiona

moreaccuratecalculationforthetwo-dimensionalvolumefraction[21]canbegivenas

80ifK>C;/<0;>>>2>111K>><2þe2D/þ2eCifCPKP0;/<0;þe12D/ifK<0;ð18ÞHð/;eÞ¼12>>2>1K>þe12D/À1ifCPKP0;/>0;>22eC>:1ifK>C;/>0

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

1e

inwhichD¼eminðjN;1Þ,C¼D2Àe2andK¼CþDje/jÀ2.TheabovesmoothedHeavisidefunctionsare2xjjNyjsymmetricto/=0returningavolumefractionof0.5.Inthispaper,foragivenfluidorasolidbody,allthecellscanbeclassifiedintothreetypes:cellswithvolumefractionlargerthan0.5arenormalcells,cellswithvolumefractionlessthan0.5butnon-zeroaresmallcells,andotherwisetheyareemptycells.Notethatanor-malcellforafluidorasolidbodycorrespondingto/>0isalsoasmalloremptycellforthefluidorsolidbodycorrespondingto/<0andviceversa.Itisnotedthatthereareotheralternativebutmoreinvolvedwaystocalculatevolumefractionsaccordingtocell-faceapertures[6].Inthispaper,theabovesimpleapproachispreferredforitssimplicityandeasyimplementation.AswillbeshowninSection7,forawiderangeoftestcasesitissufficienttoproduceverysatisfactoryresults.4.Amixingprocedure

Inthepresentmethod,theconservativequantitiesinallnormalcells,smallcellsandnewlycreatedorvan-ishedemptycellsareupdatedbyeitherEq.(9)or(8).Forasmalloremptycell,astablefluidstatemaynotbeobtainedbasedonthetimestepcalculatedaccordingtothefullgridsizeCFLcondition.Therefore,itissug-gestedthatthefluidinthosecellsshouldbemixedwiththatoftheirneighboringcells.Asthetargetedneigh-boringcellsarepreferredtobenormalcells,themixingdirectionischosenjustoppositetothatoftheextendingequation(14).Similarly,thetargetedcellsaredeterminedfromthelevelsetnormaldirection.Sup-posecell(i,j)isasmallcell,asshowninFig.3,forthefluidcorrespondingto/>0,thetargetedcellinthexdirectionischosenascell(i+1,j)ifNxi;j>0otherwisecell(iÀ1,j).Thetargetedcellintheydirectionis

(i-1,j+1)(i,j+1)Targetcellin the y directionNi,jNi,jyxSmall cell(i-1,j)Ni,j(i,j)Target cellin the x directionFig.3.Thetargetcellsforasmallcelloffluidcorrespondingto/>0.560X.Y.Huetal./JournalofComputationalPhysics219(2006)553–578

chosenascell(i,j+1)ifNyi;j>0otherwisecell(i,jÀ1).Forthefluidassociatedwith/<0,thecorrespondingoppositetargetedcellsarechosen.Thechangesoftheconservativequantitiesinthexandydirectionsforthesmallcellsandthetargetedcellsarecalculatedby

Mxi;j

¼ÀMxxtrg

ÂÃbxi;j

¼ðaxtrgUxtrgÞai;jÀðai;jUi;jÞaxtrg;ai;jþaxtrg

hibyi;j

¼ðaytrgUytrgÞai;jÀðai;jUi;jÞaytrg;ai;jþaytrg

ð19Þ

y

Myi;j¼ÀMytrg

wherethesubscriptsxtrgandytrgaretheindicesoftargetedcellsinthexandydirections,respectively.Notethatconservationismaintainedsincetheconservativequantitiesasmallcellobtainedfromatargetcellcorresponds

y

toalossofthetargetcell.Thetermsbxi;jandbi;jarethefractionsofmixingintherespectivexandydirectionsyandsatisfybxi;jþbi;j¼1.Asimpleandrobustchoiceistosetthefractionto1forthedirectionwithlargernor-ymaldirectioncomponent,thatis,maxðjNxi;jj;jNi;jjÞ.Onecanalsochooseanothersmootherbutslightlydissipa-2y2ytiveformbydefiningx-directionmixingfractionasbx¼jNxi;jjandy-directionmixingfractionwithb¼jNi;jj.Accordingtoournumericalexperiencestheformergivesbetterresultsforcomplexstaticboundaryproblemsandthelatterseemsbettersuitedformulti-fluidproblemsandcomplexboundaryproblemswithmovingbodies.Notethatotherwaystodefinemixingfractionsarepossible.Itisevenpossibletoconsiderthecorre-spondingcellinthediagonaldirection(forexample,thecell(iÀ1,j+1)inFig.3)asatargetcellbyintroducing

xxyy

anewmixingfractionbxyi;jwithrelationbi;jþbi;jþbi;j¼1andinsertinganewmixingtermintoEq.(3).

Afterthe(mixing)changesarecalculatedtheconservativequantitiesforonefluidinthenearinterfacecellsareupdatedby

󰀄󰀅ÃXX

xnþ1nþ1nþ1nþ1

þMþMy;ð20Þai;jUi;j¼ai;jUi;j

k

Ã

l

þ1nþ1

whereðani;jUi;jÞaretheconservativequantitiesattimestepn+1beforemixing.Here,thesecondandthirdtermsontheright-handsidearethesummationstakenforallthemixingoperationsinthexandydirections,respectively.Actually,theindiceskandlinEq.(20)neednotbedetermined.AstheexchangesquantitiesateverymixingoperationsofEq.(19)isaccumulatedforeverysmallandtargetcells,oneneedonlyupdatetheseaccumulatedquantitiesoncebyEq.(20)forallcellsinnearinterfaceregion.Easytofindthatonlythosecellswithnon-vanishingexchangedquantitiesareeffected.Toavoidpossibleinstability,onlyforthenormalcellsthefluidstatesareupdatedfromcell-averagedconservativequantities.Fortheothercells,thecorrespondingfluidstatesareupdatedbytheextendingequation(14)ifneeded.Notethatthepresentmixingproceduretreatsvanishedandnewlycreatedemptycellautomatically.Fortheformercase,theresidualconservativequantitiesarealltransportedtotargetcellsand,forthelattercasealltheconservativequantitiesinanewlycreatedsmallcellaretransportedfromitstargetcells.

5.Interfaceexchanges

Toobtainthemomentumandenergyexchangesacrosstheinterface,theproposedRiemannproblemsasso-ciatedwithinterfaceinteractionsaresolvedonthegridpointswithinnearinterfaceregionbandalongtheinterfacenormaldirection.WithEq.(3),fortheinteractionsbetweenfluidandacomplexboundarytheinter-facevelocityisdeterminedfromagivenprescribedorinertiallycoupledsolidbodyvelocityvI=(uI,vI)by

vI¼ðvrgÁNÞN:

ð21Þ

Forproblemswithonlyweakormoderateshocksandexpandingwaves,simplenon-iterativeapproximateRiemannsolvers[44]aresufficientforreasonableresults,otherwisemoreaccuratesolversareneeded.Inthispaper,theinterfaceconditionisobtainedbytheinterfaceinteractionmethodofHuandKhoo[20]withoutmodification.

AftertheinterfaceinteractionhasbeensolvedtheinterfacepressurepIandthenormalvelocityvI”(uI,vI)areobtained.Hence,forthefluidcorrespondingto/>0,themomentumandenergytransferredtoitare

^PðDCÞ¼pIDCNIX

and

^EðDCÞ¼pIDCNIÁvI;X

ð22Þ

X.Y.Huetal./JournalofComputationalPhysics219(2006)553–578561

(i,j+1)Ai,j+1/2i,jAi-1/2,jNi,jAi+1/2,j=0(i,j)(i+1,j)(i-1,j)NIAi,j-1/2(i,j-1)Fig.4.Schematicofconservativediscretizationforacutcell.respectively.HereDCandNIaretheinterfacesegmentlengthorareaandnormaldirection,respectively.^P¼ðX^P;X^PÞstandsforthetransferredmomentumintherespectivexandydirectionsandX^EstandsforXxy

thetransferredenergy.Inthispaper,atwo-dimensionalapproximationofEq.(22)isusedas

^PXxðDCi;jÞ¼pIðAiþ1=2;jÀAiÀ1=2;jÞDx;^PðDCi;jÞ¼pIðAi;jþ1=2ÀAi;jÀ1=2ÞDy;X

y

ð23Þ

and

^PðDCi;jÞþvIX^PðDCi;jÞ:^EðDCi;jÞ¼uIXXxy

AsshowninFig.4,thenormalizedinterfacesegmentlengthinthecellisapproximatedas

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

22

DCi;j¼ðAiþ1=2;jÀAiÀ1=2;jÞþðAi;jþ1=2ÀAi;jÀ1=2Þ:

ð24Þ

ð25Þ

ToenforceconservationtheinterfacenormaldirectionNIinEq.(22)isapproximatedwiththecell-faceaper-turedifferences.Strictly,NIhereisnotthesameasthelevelsetnormaldirectionatthegridpointNi,j(seeFig.4)whichisalsousedtoapproximatetheinterfacedirectionwhensolvingfortheinterfaceconditionforthesamecell.Thelatteriscalculatedbyastandardcentraldifferencefromthegridpointlevelsetvaluesandisalsousedtocalculatecellvolumefractions,extendingquantities,andformixingsmallcells.Thisapparentinconsistencycanbecorrectedbyreplacingthelevelsetnormaldirectionwiththeformeratthesegridpoints.However,asbothformulasareapproximationswiththesameorder,theseslightdifferenceshavealmostnoeffectonthere-sults.Notethatforthefluidcorrespondingto/<0thetransferredquantitieshaveoppositesigntothatofEqs.(23)and(24).Iftheinterfaceseparatestwofluidsthemomentumandenergytransferredbetweenthetwofluidshavethesamevaluesbutoppositesignandthereforesatisfytheoverallconservationproperty.6.Implementation

Ageneralprocedureoftheconservativeinterfacemethodcanbesummarizedasfollows:

1.Giventhecell-averageddensityofconservativequantitiesUnandfluidstatesforallgridpointswithinnor-malcells,extendthefluidstatestothegridpointscorrespondingtothesmalloremptycellsinnearinterfaceregionontheothersideoftheinterfaceviaEq.(14).

562X.Y.Huetal./JournalofComputationalPhysics219(2006)553–578

2.SolvetheinterfaceconditionsonthegridpointsinthenearinterfaceregionwiththeinterfaceinteractionmethodofHuandKhoo[20],andcalculatethemomentumandenergyexchangesacrosstheinterfacewithEqs.(23)and(24).

3.Proceedtothenexttimestepandcomputetheunmixedconservativequantities(an+1Un+1)*inallfullcellsandthecellsinnearinterfaceregionviaEq.(8).

4.Calculatethelevelsets/n+1inthenearinterfaceregionforthenewtimestepviaEq.(12)withtheproce-dureoutlinedin[20]whichupdatethelevelsetwithinterfacevelocityratherthanthefluidvelocity.

5.Re-initializethelevelsetsfortheentiredomainexceptthenearestregionoftheinterface,usuallythefirstnearestcell-layer,viaEq.(13)(detailscanbefoundin[20]).

6.From/n+1calculatethecell-faceaperturewithEq.(16),thelevel-setnormaldirectionwithstandardcentraldifferencesandvolumefractionan+1inthenearinterfaceregionbyEqs.(17)or(18).

7.InvokethemixingprocedureforthesmallcellsandnewlycreatedorvanishedemptycellswithEqs.(19)and(20)toobtainan+1Un+1from(an+1Un+1)*.

8.Updatethecell-averageddensityofconservativequantitiesUn+1andfluidstatesforallgridpointswithinnormalcells.Formulti-fluidproblems,exceptsteps4–6,allstepsarecalculatedseparatelyforeachfluid.Forcomplexstaticboundaryproblemsthesteps5and6maynotbenecessary.Ifthesolidbodyassociatedwiththecomplexgeometrymoveswithaprescribedvelocityvrg(t),thelevelsetmoveswiththisprescribedsolidbodyvelocity.Inthecasethatthesolidbodymovesbyinertialcouplingthesolidbodyvelocityisupdatedfromthetotalmomentumvariation.AccordingtoEq.(4),theaverageaccelerationinonetimestepis

1X^P

XðDCi;jÞð26Þarg¼

mrgi;jandthevelocityofthesolidbodyafteronetimestepcanbeobtainedfrom

þ1nvnrg¼vrgþargDt:

ð27Þ

ForhigherordertimeaccuracyaRunge–Kuttatimeintegrationcanbeused.IneachRunge–Kuttasub-stepalltheabovestepsareinvokedonce,exceptforthere-initializationstep5whichisdoneonlyonceafterthelastsub-step.Currently,onlytwo-dimensionalproblemsareconsidered.However,theextensionofthepresentmethodtothreedimensionsisfairlystraightforward.TheonlymoreevolvedworkisthecalculationofcellfaceaperturesandtheapproximateformoftheinterfacesurfacepatchDCforathree-dimensionalcutcell.7.Numericalexamples

Thefollowingnumericalexamplesareprovidedtoillustratethepotentialofthepresentinterfacemethod.Foralltestcases,one-phasecalculationsarecarriedoutwithafifth-orderWENO-LF[24]andathird-orderTVDRunge–Kuttascheme[41].Forone-dimensionalexamples,thenumberofgridpointsis200andtheref-erencedexactsolution,ifexists,issampledon1000gridpoints.Inthetwo-dimensionalexamples,ifnotmen-tionedotherwise,theupperandlowerboundariesareimposedwithreflecting-wallboundaryconditionandtheleftandtherightboundariesimposedwithoutflowboundaryconditionwithconstantextrapolation.AllthecomputationsarecarriedoutwiththeCFLnumberof0.6.7.1.Gas–gasinteraction(I)

CaseI-A:Weconsideragas–gasinteractionproblemwiththefollowinginitialconditions

&

ð1;0:5sin½pð1À0:5xÞ󰀄;1;1:4Þifx<0:5;

ðq;u;p;cÞ¼

ð1;0:5sin½pð1À0:5xÞ󰀄;1;1:8Þifx>0:5

ð28Þ

andreflectingwallboundaryconditionappliedatbothx=0andx=1.Thecaseiscomputeduptotimet=1.Weexaminethenumericalsolutionswiththereference‘‘exact’’solutioncomputedwith1600gridpoints.Fig.5showsthecalculatedvelocityanddensityprofilesattimet=0,0.25,0.75and1.Attheearlystageof

X.Y.Huetal./JournalofComputationalPhysics219(2006)553–578

22563

t=01.751.51.25ExactdensityExactvelocityCalculateddensityCalculatedvelocityt=0.251.751.51.25ExactdensityExactvelocityCalculateddensityCalculatedvelocityDensity & Velocity10.750.50.250Density & Velocity00.250.50.75110.750.50.250-0.25-0.5-0.25-0.5(a)21.751.51.25x(b)200.250.5x0.751t=0.75ExactdensityExactvelocityCalculateddensityCalculatedvelocityt=11.751.51.25ExactdensityExactvelocityCalculateddensityCalculatedvelocityDensity & Velocity10.750.50.250Density & Velocity00.250.50.75110.750.50.250-0.25-0.5-0.25-0.5(c)x(d)00.250.5x0.751Fig.5.Gas–gasinteraction:CaseI-A.interactiontheflow,asshowninFigs.5bandc,iscontinuouswithprimarilyonlyexpansionandcompressionwavestransmittedacrosstheinterface.Latertheflowdepictsthepresenceofshockdiscontinuities,afterthetwocompressionwaveshavebeenreflectedfrombothendsandinteractwiththeinterfaceandwitheachother(seeFigs.5candd).Theseshockwavespassthroughtheinterfaceataboutt=0.and0.80,respectively.Onecanobservethatthecomputedresultsareingoodagreementwiththe‘‘exact’’solution.Resolutionstudiesarecarriedouttomeasurethenumericalconvergencerateofthemethod.Fordifferentresolutions,Fig.6givesthePlllvariationsoftotalentropyoftheleftmediumsl¼ieliai,whereeiandaPiaretheentropyandthevolumefractionoftheleftmediumincelli,respectively,andoftotalentropysT¼ieiai,whereeiandaiaretheen-tropyandvolumefractionoftheleftorrightmediumincelli,respectively.WemeasuretheentropyerrorsEs=|stÀsexact|/sexactatt=0.2,0.and0.8correspondingtothetimewhentheentropyassumestheextremevalueinacontinuousfloworwhentheshockwavespasstheinterface.TheentropyerrorsElsandETsfortheleftandtotalmedium,respectively,andtheaveragedorderofconvergenceRcaregiveninTable1.Onecanfindthathigh-orderaccurateresultsareobtainedforthecontinuousflow.However,whenthereisashockwaveandespeciallywhentheshockwavepassestheinterface,largererrorsareproducedwithanorderofcon-vergencedecreasingtowardsunity.

CaseI-B:Weconsiderawellknownair–heliumshocktubeproblemwiththefollowinginitialconditions

&

ð1;0;1;1:4Þifx<0:5;

ðq;u;p;cÞ¼ð29Þ

ð0:125;0;0:1;1:667Þifx>0:5:

5X.Y.Huetal./JournalofComputationalPhysics219(2006)553–578

0.020.080.070.06Exact50gridpoints100gridpoints200gridpoints400gridpoints0.0180.0160.014Exact50gridpoints100gridpoints200gridpoints400gridpointsEntropy (leftmedium)0.050.040.030.020.010-0.01-0.02-0.0300.250.50.751Entropy (total)0.0120.010.0080.0060.0040.002000.250.50.751(a)time(b)timeFig.6.CaseI-A:convergencetest.Table1

ErrorsandconvergencerateforCaseI-A1/h50100200400Rc0:2Ets¼À

¼0:2

EtsT

0:Ets¼À

¼0:

EtsT

0:8Ets¼À

¼0:8

EtsT

3.7·10À42.6·10À54.9·10À.0·10À73.3

4.1·10À49.2·10À52.4·10À55.7·10À62.5

5.2·10À32.1·10À37.1·10À41.2·10À41.8

4.9·10À32.1·10À39.0·10À43.5·10À41.3

6.7·10À33.1·10À31.1·10À32.8·10À41.6

5.6·10À32.5·10À31.1·10À34.9·10À41.2

Typicalcomputationalresultsattimet=0.15areshowninFig.7.Onecanfindthecalculatedresultsareinreasonablygoodagreementwiththeexactsolution.Inparticular,theyhavealmostidenticalandcorrectshockstrengthandspeed.Theinterfacepositionisalsocapturedaccurately.Notethattherearelargerdensitydis-crepanciesneartheinterface(alsoseeninthefollowingCaseI-C)comparedwithtotheresultsofFedkiwetal.[12](seetheirFig.8)andofHuandKhoo[20](seetheirFig.3).Thisisbecausetheisobaricfixusedintheirworksproducesbetterdensityprofilesnearinterfaces.Similartreatmentscouldalsobeemployedinourmeth-odbymodifyingthestatesbeforecomputingtheinterfaceconditions.However,weprefernottousethesaidfixsinceitmayviolateconservationconsiderably.AsshowninFig.7d,thecalculatedinterfacelocationsug-gestsbetterconvergencepropertythanthatofFedkiwetal.[12](seetheirFigs.7and8),thelattershowsthattheinterfacepositionappearstobeoffbyonecellforalllevelsofrefinementasconservationattheinterfaceisnotsatisfied.

Wecomputetherelativemassandenergyvariationsforeachfluidseparatelyduringthecomputationby

Pi¼Nnn

0aUiDx;ð30ÞV¼Pii¼¼N00aUDxi¼0iwhereUiisthecell-averagedmassdensityorenergydensityoncelli.Thesuperscriptsnand0representthenth

timestepandtheinitialcondition,respectively.Inasimilarway,thetotalmassandenergyvariationsaregivenas

hPiPi¼Nnni¼Nnn

ði¼0aUiÞleftþði¼0aUiÞrightDx

i:ð31ÞVtotal¼hPPi¼N00i¼N00

ði¼0aUiÞleftþði¼0aUiÞrightDxFig.8comparesmassandenergyvariationofeachfluidtothoseobtainedbytheinterfaceinteractionmethodofHuandKhoo[20](denotedasI-GFM).Noconservationerrorisobservedforthecurrentmethodastheerrorisoftheorderofmachineaccuracyof10À14.

X.Y.Huetal./JournalofComputationalPhysics219(2006)553–578

1ExactCalculated565

1ExactCalculated0.80.8Pressure0.6Density0.250.50.7510.60.40.40.201ExactCalculated0.200.250.50.751(a)x(b)x1.80.8ExactCalculated:200pointsCalculated:400pointsGammaVelocity0.61.60.41.40.20(c)00.250.5x0.7511.20.6250.630.635(d)x0.0.5Fig.7.CaseI-B:air–heliumshocktubeproblem.CaseI-C:Wecomputeastiffair–heliumshocktubeproblem,takenfrom[1,20],wheretheinitialpressuredifferenceismuchlargerthanthatinCaseI-B.Theinitialdataare

&

ð1;0;500;1:4Þifx<0:5;

ðq;u;p;cÞ¼ð32Þ

ð1;0;0:2;1:667Þifx>0:5:Thecomputedresultsattimet=0.015areshowninFig.9.Thereisgoodagreementwiththeexactsolution.Attherightendoftherarefactionwavesonthepressureandvelocityprofiles,onecanobservemildover-shootswhichbecomessignificantin[1](seetheirFig.6).Theirschemealsoproducesmorenumericalviscositywhichleadstostrongersmearingoftheshockfront.Itrequiresamuchfinerdistributionofabout800gridpointsforacomparablysharpshockfront.NotethattheseovershootsarereplacedbyundershootsintheresultsofHuandKhoo[20](seetheirFig.4).

CaseI-D:2Dair–heliuminteraction.Inthis2-Dproblem,wecomputeaMach6airshockwaveinteractingwithacylindricalheliumbubble.Numericalcomputationsforthesameproblemcanbefoundin[3]whotreatairandheliumasthesamefluidtoavoiddifficultiesattheinterface.Forthiscase,theinitialconditionsare

8

ðq¼1;u¼0;v¼0;p¼1;c¼1:4Þpre-shockedair;>>>><ðq¼5:268;u¼5:752;v¼0;p¼41:83;c¼1:4Þpost-shockedair;

ð33Þðq¼0:138;u¼0;v¼0;p¼1;c¼1:667Þheliumbubble;>>qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi>>:2

/¼À0:025þðxÀ0:15Þþy2levelset;where/60representstheheliumand/>0representstheair,depictingaheliumbubbleofradius0.025at

(0,0.15)whichistobeimpactedbyashockwaveinitiatedatx=0.1.ThecomputationhasbeencarriedoutwiththreeincreasingresolutionsofDx=Dy=5·10À3,2.5·10À3and1.25·10À3.

566X.Y.Huetal./JournalofComputationalPhysics219(2006)553–578

1.0081.002I-GFMThispaperI-GFMThispaper1.006Mass (leftmedium)1.001Mass (rightmedium)0.050.10.151.00411.0020.99910.998000.050.10.15timetime1.003I-GFMThispaper1.003I-GFMThispaper1.0021.0021.001Energy (total)0.050.10.15Mass (total)1.001110.9990time0.99900.05time0.10.15Fig.8.CaseI-B:conservationtest.Fig.10showsthecalculateddensitycontoursatfourselectedtimeinstances(Dx=Dy=1.25·10À3).TheseresultsshowafairlygoodagreementwiththoseofBagabirandDrikakis(theirFig.6).Thesecondaryreflectedshockwaveandtriple-waveconfigurationsarecalculatedwithgoodresolutionasshowninFigs.10b–d.AsBagabirandDrikakis’scomputationalmodeldoesnottreatinterfaces,theirbubbleandjetshapesaresmearedandhardtoidentifyespeciallyafterthetransmittedshockwavehaspassedthebubble.Inthepresentcalcu-lationtheseshapechangesareaccuratelycaptured,includingthestronginstabilitiesattheinterfacewhichevenproduceisolatedairdropletsintheheliumbubble(seeFigs.10bandc).Alsowellrepresentedistheinter-faceafterthejetcollideswiththedownstreambubblesurfacewhichevolvestoproducesmallseparatedheliumpocketsinthebackgroundairmedium(seeFig.10d).Anotherfeatureofthepresentmethodisthatthecom-plextopologicalchangesistreatedrathernaturallywithoutdifficulty.TheinterfaceatthreeselectedtimeinstancesisshowninFig.11a.Itisobservedthattheheliumbubbleevolvesintoverycomplexgeometrieswithmanyseparatedparts.Inthepresentmodel,whentheseparateddroplet/pocketbecomessmallerthanthegridsize,theunresolvedmaterialisautomaticallydiscardedasthere-initializationprocedureignoresdetailssmal-lerthanthegridsizeandreconstructsthenewlevelsetsfromtheresolvedregion.ThesearefurthershowninFig.11bwhichsuggeststhedeletingofunresolvedheliumbyaseriesofstair-step-likechangesoftotalheliummass.Notethatthesemasschangeshappenlongtimeafterthefirsttopologicalchanges(ataboutt=1.1·10À2)suggestingthattopologicalchangesdonotviolatetheconservationpropertyofthecurrentmethod.Fig.11cshowsthetotalforceactingontheupper-halfoftheheliumbubble.Itisobservedthatsec-ondorderorsuper-linearconvergenceforthepeakforcesbeforethefirsttopologicalchangeatabout

X.Y.Huetal./JournalofComputationalPhysics219(2006)553–578

500ExactCalculated567

4.543.53ExactCalculated400PressureDensity0.250.50.7513002.521.520010010.502ExactCalculatedExactCalculated(a)015x(b)0.250.5x0.751121.8VelocityGamma00.250.50.75191.661.431.201(c)x(d)00.250.5x0.751Fig.9.CaseI-C:stiffair–heliumshocktubeproblem.t=1.1·10À2.Afterthetopologicalchanges,itishardtoevaluatetheconvergencerateastheinterface-linelengthisstronglydependentoninterfaceinstabilities,whichsuggestnosimpledependencyongridsizes.Fig.11dcomparesthejetandupstreambubblesurfacetrajectoriescalculatedwithdifferentgridsizes.Asec-ond-orderconvergencerateisobtainedbycomputingtotalrelativeerrorsbetweensuccessivegridsizes.7.2.Gas–waterinteraction(II)

CaseII-A:Gas–watershocktubeproblem.Thisproblemistakenfrom[20].ArealgasequationofstateisusedforthegasandTait’sequationofstateforthewatermedium.Theinitialconditionsaregivenas

&

ð0:01;0;1000;2Þifx<0:5;

ðq;u;p;cÞ¼ð34Þ

ð1;0;1;7:15Þifx>0:5:Inthisproblemthehighpressuregasexpandsslowlycomparingtothetransmittedandreflectedwavefront

speeds.Fig.12showsthecomputedresultsattimet=0.0008.NotethatthepresentresultsindicateslightlymoreaccuratesolutionsfortherarefactionwaveingasmediumandthetransmittedshockwaveinwaterthanthoseofHuandKhoo[20];thelatterproducesaslightovershotattherarefactionwaveandsmalldifferencesinthevelocityatthewater–gasinterface(seetheirFig.12).

CaseII-B:Shockimpactingonanair–waterinterface.Wehereconsidertheinitialconditions

&

ð1:037578;6:0151;1000;7:15Þifx<0:7;

ðq;u;p;cÞ¼ð35Þ

ð0:001;0;1;1:4Þifx>0:7:Thisisanairbubblecollapseprobleminonedimension.Whilethegasbubblecollapses,theunderwater

shockwaveimpactsattheinterfaceresultinginaweakshockwavetransmittedintotheairandaverystrong

568X.Y.Huetal./JournalofComputationalPhysics219(2006)553–578

0.040.030.020.010-0.01-0.02-0.03-0.04t=6x10-30.040.030.020.010-0.01-0.02-0.03-0.04t=1.2x10-2y(a)0.10.1250.15x0.1750.2y(b)0.10.1250.15x0.1750.20.040.030.020.010-0.01-0.02-0.03-0.04t=1.8x10-20.040.030.020.010-0.01-0.02-0.03-0.04t=2.4x10-2y(c)0.1250.150.175x0.20.225y(d)0.1750.20.225x0.250.275Fig.10.CaseI-D:35contoursofdensity,from0to9forselectedtimeinstances.rarefactionwavereflectedbacktothewatermedium.Thecalculatedresultsattimet=0.003areshowninFig.13.Onecanfindthat,whileinoverallagreementwiththeexactsolution,theresultsproducedwiththepresentmethodpredictamoreaccuratepressureprofilethanthatcalculatedwiththeinterfaceinteractionmethodofHuandKhoo[20](denotedasI-GFMinFig.13).NotethattheI-GFMisalsoafairlyrobustmethod.AsthepressureprofileinFig.13aisplottedinlogarithmscale,theoscillationproducedbyI-GFMisstillverylimitedinmagnitudeanddoesnotcauseanystabilityproblem.OurnumericalexperimentssuggestthatI-GFMstillabletoproducereasonableresultsalthoughnotbeingexactlyconservative.

CaseII-C:Underwaterexplosion.Weconsideranunderwaterexplosionproblem[45]inwhichtheexplosionproductsaremodeledasahighpressureunderwatergasbubble.Thegrowthandcollapseofanunderwatergasbubblecanbedescribedbytheone-dimensionalsphericallysymmetricconservationlaws[26]inwhichthesphericaleffectsaretakenintoaccountassourceterms.ThereareanumberofprevioussimulationsofthisproblemviaarbitraryLagrangian–Eulerian(ALE)methods[29,37,42]oranEulerianmethodwithasimplifiedmodelneglectingthedynamicsofgasbubbleandshockinterfaceinteractions[25].Thisproblemisabigchal-lengeforanyEulerianmulti-fluidmethodsincethereareverylargegasbubblevolumechangesataverylowinterfacespeedswhichheretobecomputedforaverylongphysicaltime.Therefore,averygoodconservationpropertyiscrucialtopredictthisphenomenoncorrectly.Theinitialconditionsforthisproblemare

&

ð1:63;0;83810;1:4Þif0ðq;u;p;cÞ¼ð36Þ

ð1:025;0;10;5:5Þif400>x>0:16

X.Y.Huetal./JournalofComputationalPhysics219(2006)553–578

6.000.04569

Differencetoinitialtotalheliummass(%)0.030.020.0105.004.00Deletingofunresolvedmass3.00yt=1.8x10-2t=2.4x10-2-0.01-0.02-0.03-0.04t=3x10-22.001.000.000.20.2250.25(a)0.3x0.2750.300.51(b)1.6Time(10)1.5-222.530.2Fx-325x90Fy-325x90Fx-650x180Fy-650x180Fx-1300x360Fy-1300x3601.4FyTime(10)-21.2jet-325x90downstream-325x90jet-650x180downstream-650x180jet-1300x360downstream-1300x3600.1FxandFy1.0jetdownstream0.00.8-0.10.6-0.2Fx-0.300.20.40.60.81-20.41.21.41.60.20.120.140.160.180.20.22(c)Time(10)(d)xFig.11.CaseI-D:(a)interfacesatthreeselectedtimeinstances,(b)timevariationofthedifferencebetweentotalheliummasstoitsinitialvalue,(c)forcesontheupper-halfoftheheliumbubblewithdifferentgridsizes,(d)jetanddownstreambubblesurfacetrajectorieswithdifferentgridsizes.andwhicharetakenfrom[37].Thereflectingwallboundaryconditionisusedfortheleftboundary.Wedonotimplementaboundaryconditiontotherightboundaryasthedisturbancedoesnotarrivethereinthecalcu-lations.Asin[37]astiffgasequationofstateisusedforwater,i.e.(cÀ1)qe=p+cp0wherec=5.5andp0=4921.15.ToobtainameshconvergedbubbleoscillationperiodthecalculationisperformedwithgridsizesofDx=0.2,0.1,0.05,0.025and0.0125,correspondingto0.8,1.6,3.2,6.4and12.8gridpoints(compu-tationalcells)fortheinitialgasbubbleof16cmradius.

Fig.14showsthecalculatedbubbleradiusandinterfacepressurehistorywithdifferentgridsizes,bothofwhichsuggestasecond-orderconvergencerate.Theconvergedresults,suchasabubbleoscillationperiodofabout0.2sandamaximumbubbleradiusatabout3mattimet=0.1ms,areingoodagreementwithpre-viousALEcalculation[37].Ourresultsalsoshowpressureoscillationsattheinterfacewhichhasbeenfoundinpreviousworks[29,30,45].NotethatthecurrentcalculationscapturecorrectphysicswithmuchfewergridpointforthegasbubblethanpreviousALEsimulations.Furthercalculations(notshownhere)suggestthatwhenTait’sequationisusedforthewatermedium,almostidenticalresultsareobtainedwithaslightlylargeroscillationperiod.Inthepresentcomputationstheso-called‘‘jaggedness’’(oroscillatorybehaviorwithlimitedmagnitude)phenomenawithcomplexsecondary,tertiaryandmoresuccessiveshockwavesproducedbyinter-faceinteractionscanbecapturedclearlywhenthegasbubbleisdescribedwithsimilarlyhighnumerical

570

1200X.Y.Huetal./JournalofComputationalPhysics219(2006)553–578

1.2ExactCalculated1.110.91000ExactCalculated8000.8PressureDensity00.250.50.7510.70.60.50.40.36004002000.20.10(a)x(b)1000.250.5x0.7516ExactCalculatedExactCalculated8Velocity3Gamma20000.250.5(c)x0.751(d)00.250.5x0.751Fig.12.CaseII-A:gas–watershocktubeproblem.resolutionaspreviousALEsimulations.Fig.15depictsthe‘‘jaggedness’’phenomenaattheearlystageoftheinteractionwhenapproximately600gridpointsareusedforthegasbubble.Again,theseresultsareinverygoodagreementswithpreviousALEsimulations[29,37,42].

CaseII-D:Collapseof2Dairbubblesinwater.Weconsidertheinteractionsbetweentwoneighboring2Dbubblesandastrongshockwave.ThisproblemhasbeenpreviouslysimulatedbyHankin[17]:two,6and3mmdiametercylindricalairbubblesinwaterareimpactedbya1.6GPashockwave.AccordingtoHankin[17],theschematicoftheproblemisgiveninFig.16awhichshowsa15·15mmcomputationaldomain.Theinitialdataare

8

ðq¼1;u¼0;v¼0;p¼1:013;c¼7:15Þpre-shockedwater;>>>>>ðq¼1:226;u¼54:28;v¼0;p¼16000;c¼7:15Þpost-shockedwater;>>><ðq¼1:2Â10À3;u¼0;v¼0;p¼1;c¼1:4Þairbubble;

ð37Þqffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

>22>>/¼À3þðxÀ7:5ÞþðyÀ8:5Þlevelsetfory>Àxþ11:5;>>>qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi>>:

levelsetfortherestregion:/¼À1:5þðxÀ3:75Þ2þðyÀ5Þ2Here/60representstheairand/>0representsthewater.Initially,theshockwavejustreachesthesmallair

bubble.150·150gridpoints,approximatelythesameasin[17],areusedforthecomputations.Figs.16b–dshowthetypicaldensitycontoursatt=0.18(1.8ls),t=0.34(3.4ls)andt=0.46(4.6ls)whichdepictgoodagreementwiththoseofHankin[17](theirFig.5).Onecanfindthat,althoughevenslightlylessgridpointshavebeenusedinthepresentcalculationtheshapesofthebubblesarecapturedinmuchmoredetail.Notethatwhenthesmallbubbleissplitbyjetsandcompressedtoascalesmallerthanthegridsizeitvanishesandaffectstherestflowfieldnomore.Furthermore,afterthesmallbubblehascollapsedandvanishedthe

X.Y.Huetal./JournalofComputationalPhysics219(2006)553–578

1.2103571

11020.8PressureDensity1010.60.4100ExactI-GFMThispaper0.2ExactI-GFMThispaper10-1(a)1200.250.5x0.751(b)000.250.5x0.75110Velocity6Gamma3ExactI-GFMThispaper2ExactI-GFMThispaper0000.250.5(c)x0.751(d)00.250.5x0.751Fig.13.CaseII-B:shockimpactingair–waterinterface.4.543.532.521.510.50Gridsize:20cmGridsize:10cmGridsize:5cmGridsize:2.5cmGridsize:1.25cm104Gridsize:20cmGridsize:10cmGridsize:5cmGridsize:2.5cmGridsize:1.25cmInterfacepressure(atm)00.511.522.5103Interfaceposition(m)10210110010-1(a)Time(10-1s)(b)00.51Time(10-1s)1.522.5Fig.14.CaseII-C:bubbleradiustimeandinterfacepressurehistoryfordifferentgridsizes.explosionwavesoriginatingfromthesecollapsepoints(seeFig.16c)leadtotwocavitationregions.Sincethepresentmodeldoesnotincludeanycavitationmodelthecomputationisstoppedwhenintheseregionseven-tuallyanegativesoundspeedarises.

572

106X.Y.Huetal./JournalofComputationalPhysics219(2006)553–578

104t=0.28mst=0.44mst=0.592mst=0.8ms105t=0.026mst=0.06mst=0.104ms103Pressure103Pressure102104t=0.234ms102101(a)1.81.600.250.50.751101x(b)1.41.200.250.50.7511.251.51.752xt=0.026mst=0.28mst=0.44mst=0.592ms1.41.2t=0.8mst=0.06mst=0.104mst=0.234ms1Density10.80.60.40.2000.250.50.751Density0.80.60.40.20(c)x(d)00.250.50.7511.251.51.752xFig.15.CaseII-C:thepressureanddensityprofilesforselectedearlytimeinstances.7.3.Complexboundaryproblems(III)

CaseIII-A:Movingwallproblem.Weconsideragasconfinedbetweentworeflectingwallsatxl=0andxr=0.5+urtwithconstantur=0.5.Theinitialconditionsare

½qðxÞ;uðxÞ;pðxÞ󰀄¼½1þ0:2cosðpÀ2pxÞ;2urð1ÀxÞ;qðxÞ󰀄

1:4

ð38Þ

inwhichtheentropys(x)=1.Thisproblemistakenfrom[13].Weexaminethenumericalsolutionsandcom-paringtothereference‘‘exact’’solutioncomputedwith1600gridpoints.Wecalculatethiscasetotimet=0.6.Fig.17presentstheresultsattimet=0,0.2,0.4and0.6.Onecanobservethatthecomputedresultsareinquitegoodagreementwiththe‘‘exact’’solution.Notethatthegridpointsoccupiedbythewallarefilledwithzeros.Resolutionstudiesarealsocarriedouttomeasurethenumericalconvergencerateofthemethod.SincetheanalyticsolutionissmooththeentropystaysPfinalP0constantthroughout.Weexaminethetotalentropyerrorfinal

errttlcalculatedbyerrttl¼jjsjÀ1jaj=jajandtheentropyerroronthemovingboundarygivenby

À1j,respectively,atthefinaltime.Table2showstheanerroranalysisoftheaboveerrors,whicherrbry¼jsfinalI

givesasecond-orderconvergencerateforthetotalerrorandaslightlybetterthanlinearconvergenceratefortheboundaryerror.TheseresultsareingoodagreementwithForrerandBerger[13]andArientietal.[2]usingnon-conservativemethods.

CaseIII-B:Shockdiffractiononanairfoil.WeconsideraMach1.5shockwavediffractiononaNACA0018airfoilwith+30°angleofattack.NumericalcomputationsforthesimilarproblemwithCartesiangridcanbefoundin[47]withthegridintersectionmethodandin[39]withanimmersedboundarymethod.Thepresentset-upoftheproblemisthesameasin[47]wherea1.8·2computationaldomainisemployed.Theinitialdataare

&

ðq¼1:4;u¼0;v¼0;p¼1:0;c¼1:4Þpre-shockedair;

ð39Þ

ðq¼2:607;u¼0:694;v¼0;p¼2:4583;c¼1:4Þpost-shockedair:

X.Y.Huetal./JournalofComputationalPhysics219(2006)553–578

1515573

t=01212t=0.019y6y63300369121500(a)15x36(b)15x91215t=0.0341212t=0.04699y6y63300036912150(c)x(d)36x91215Fig.16.CaseII-D:50contoursofdensityfrom0to1.3atdifferenttimeinstances.Thelevelsetisbuiltfromapolylineconsistingof200controlpoints[21].Fig.18ashowsthecomputedresultsona360·400gridattimet=0.whichsuggestsagoodagreementwiththeshockstructurefromtheexper-imentalresultsofMandellaandBershader[31]atanapproximatelysimilartime.ComparedwiththeresultsofRipleyetal.[39]obtainedwitha500·500grid,thepresentmethodshowshigherorderaccuracybywhichamuchsharpershockwavefrontandrefinedflowstructuresneartheairfoilheadandtipareproduced.Notethatthediscontinuitiesnearthesolidbodyaresharplypredictedwhicheliminatestheproblemsassociatedwithartificialnumericalboundarylayerfoundin[38,47].Byimposingthefreestreamvelocityanddensitywith

F

thepost-shockvalueswecalculatetheliftanddragcoefficientsbyCl¼1y2,Cd¼1Fx2whereFxandFyarethetotalmomentumexchangerateintegratedalongtheboundaryinthexandydirections,respectively,andListheairfoilchordlength.Thecalculatedliftanddragcoefficientswithdifferentgridsizesareshownin

2q1u1L

2q1u1L

574

1.25X.Y.Huetal./JournalofComputationalPhysics219(2006)553–578

t=0ExactdensityExactvelocityCalculateddensityCalculatedvelocity1.25t=0.2ExactdensityExactvelocityCalculateddensityCalculatedvelocity11Density & Velocity0.75Density & Velocity00.250.50.7510.750.50.50.250.250000.250.50.751(a)1x(b)xt=0.4ExactdensityExactvelocityCalculateddensityCalculatedvelocityt=0.60.75ExactdensityExactvelocityCalculateddensityCalculatedvelocity0.75Density & VelocityDensity & Velocity0.50.50.250.25000.250.50.7510(c)x(d)00.250.5x0.751Fig.17.CaseIII-A:movingwallproblem.Table2

ErrorsandconvergencerateforCaseIII-A1/herrttlerrbry502.9·10À44.06·10À31006.87·10À51.74·10À32001.71·10À57.7·10À44004.31·10À63.4·10À48001.07·10À61.5·10À4Rc2.031.2

Fig.18bwhichsuggestaboutsecond-orderconvergencerateforthepeakliftcoefficientandlinearconvergenceforthedragcoefficient.

CaseIII-C:Theliftoffofacylinderbyashockwave.Thisexampleofaflowwithcomplexmovingboundaryistakenfrom[11]ontheliftoffofacylinderbyashockwave.Arigidlight-weightcylinder,initiallyrestingonthelowerwallofatwo-dimensionalchannel,isdrivenandliftedupwardsbyaMach3shockwave.Theprob-lemhasalsobeensimulatedwithnon-conservativemethodsbyForrerandBerger[13]andbyArientietal.[2]whichimplementthemovingboundaryconditionwithmirrororinjectionoperations.Theinitialdataforthisproblemare

8

ðq¼1;u¼0;v¼0;p¼1Þpre-shockedair;>><

ðq¼3:857;u¼2:629;v¼0;p¼10:333Þpost-shockedair;ð40Þqffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi>>:/¼À0:05þðxÀ0:15Þ2þðyÀ0:05Þ2levelset:

X.Y.Huetal./JournalofComputationalPhysics219(2006)553–578

21.7521.51.2512.5Cd-90x100Cl-90x100Cd-180x200Cl-180x200Cd-360x400Cl-360x400575

CdClandCd(10)-21.5y0.750.51Cl0.50.2500(a)00.5x11.5(b)0123time4567Fig.18.CaseIII-B:(a)35fluidpressurecontoursfrom0to3.5,(b)timedependentliftanddragcoefficients.0.2t=00.150.1y0.050(a)0.200.250.50.75xt=0.140.150.1y0.050(b)0.200.250.50.75xt=0.2550.150.1y0.050(c)00.250.50.75xFig.19.CaseIII-A:60contoursoffluidpressurefrom0to28forselectedtimeinstances.576X.Y.Huetal./JournalofComputationalPhysics219(2006)553–578

Here,/60representsthecylinderwitharadiusof0.05anddensityof7.6and/>0representstheairinthechannel.Thecomputationaldomainis1.0·0.15,forwhichtheleftinletboundaryassumesthepost-shockparameters.DifferentgridsizesofDx=Dy=2.5·10À3,1.25·10À3and6.25·10À4areused.Fig.19showsthepressurecontoursatthreedifferenttimest=0.0,0.014and0.0255,whichcorrespondtotimest=0,500and910ofFalcovitzetal.[11],respectively.Ourresultsshowthatattimet=0.028(t=1000in[11])thecylinderhasalreadyleftthetopboundaryofthedomain.TheseresultsagreewithFor-rerandBerger[13]andArientietal.[2].AstrongvortexappearsunderthecylinderintheresultsofForrerandBerger[13](seetheirFig.5).However,asshowninFigs.19bandc,thisvortexcannotbeenfoundhere.Thesamephenomenonisfoundafterwehaveexaminedourlowresolutionresults(notshownhere)andtheresultsin[2](theirFig.19).Onepossiblereasonforthisdiscrepancymaybecausedbythespace–timesplit-tingschemeemployedin[13].Thedifferenttreatmentofthetangentialdirectionvelocitiesfortheghostnodesmayeffectsthenumericaldissipationconsiderably.Inthepresentresults,themeasuredcylindercen-tersatt=0.0255withincreasingresolutionare(0.659,0.132),(0.9,0.145)and(0.1,0.147),correspondinggridsizesareDx=Dy=2.5·10À3,1.25·10À3and6.25·10À4.Thereresultsimplyasuper-linearconver-gencerate,slightlybetterthanthatofArientietal.[2].Thisisperhapstobeexpectedsincethepresentcon-servativeinterfacemethodmayimprovethesolutionespeciallywhenthecalculationiscarriedoutonlongtimescales.

8.Concludingremarks

Inthispaper,wehavedevelopedaconservativeinterfacemethodsuitableforbothmulti-fluidproblemsandcomplexboundaryproblems.Theconservationpropertyremovesorreducesanyconservativeassociatederrorsespeciallywhentherearestrongorlongphysicaltimescaleinteractionsacrosstheinterface.AsthemethodisconstructedbasedonastandardCartesianfinitevolumemethodandlevelsettechnique,itmain-tainsthesimplicityofGFM-likemethodsforimplementationandhandlestopologicalchangesnaturally.Thepresentmethodalsooffersafairlysimplewayofimplementationinmulti-dimensionandformulti-leveltimeintegrationswithoutspace–timesplitting.Anumberofnumericalexamplesinonedimensionarestudiedwithcomparisonstoexactsolutionswhiletwo-dimensionalproblemsarecalculatedandcomparedtoexperimentsandresultsofpreviousworks.Theobtainedresultssuggestthatthepresentmethodexhibitsalargedegreeofrobustnessandaccuracywithquitegoodconvergenceproperties.Asthedifferenttypesofinterfaceproblemsinthisworkaretreatedwiththesameapproachourmethodmayprovideawaytowardsuniversalmethodsforinterfaceproblems.References

[1]R.Abgrall,S.Karni,Computationsofcompressiblemultifluids,J.Comput.Phys.169(2001)594–623.

[2]R.Arienti,P.Hung,E.Morano,J.E.Shepherd,AlevelsetapproachtoEulerian–Lagraniancoupling,J.Comput.Phys.185(2003)213.

[3]A.Bagabir,D.Drikakis,Machnumbereffectsonshock–bubbleinteraction,ShockWaves11(2001)209–218.

[4]J.A.Benek,T.L.Donegan,J.L.Steger,Extendedchimeragridembeddingschemewithapplicationtoviscousflows,AIAAPaperNo.87-1126,1987.

[5]M.J.Berger,R.J.LeVque,AnadaptiveCartesianmeshalgorithmfortheEulerequationsinarbitrarygeometries,AIAAPaperNo.-1930,19.

[6]R.Caiden,R.P.Fedkiw,C.Anderson,Anumericalmethodfortwo-phaseflowconsistingofseparatecompressibleandincompressibleregions,J.Comput.Phys.166(2001)1–27.

[7]Y.L.Chiang,B.vanLeer,K.G.Powell,SimulationofunsteadyinviscidflowonanadaptivelyrefinedCartesiangrid,AIAAPaperNo.92-0443,1992.

[8]M.H.Chung,AlevelsetapproachforcomputingsolutionstoinviscidcompressibleflowwithmovingsolidboundaryusingfixedCartesiangrids,Int.J.Numer.Meth.Fluids33(2001)1121–1151.

[9]P.Colella,Multidimensionalupwindmethodforhyperbolicconservationlaws,J.Comput.Phys.87(1990)171–200.

[10]S.F.Davis,Aninterfacetrackingmethodforhyperbolicsystemsofconservationlaws,Appl.Numer.Math.10(1992)447–

472.

[11]J.Falcovitz,G.Alfandary,G.Hanoch,Atwo-dimensionalconservationlawsschemesforcompressibleflowswithmoving

boundaries,J.Comput.Phys.72(1997)449–466.

X.Y.Huetal./JournalofComputationalPhysics219(2006)553–578577

[12]R.Fedkiw,T.Aslam,B.Merriman,S.Osher,Anon-oscillatoryEulerianapproachtointerfacesinmultimaterialflows(theghostfluid

method),J.Comput.Phys.152(1999)457–492.

[13]H.Forrer,M.Berger,FlowsimulationonCartesiangridsinvolvingcomplexmovinggeometriesflowsInt.Ser.Numer.Math.,vol.

129,Birkha¨user,Basel,1998.

[14]H.Forrer,Jeltsch,Ahigh-orderboundarytreatmentforCartesian-gridmethod,J.Comput.Phys.140(1998)259–277.

[15]J.Glimm,D.Marchesin,O.McBryan,Anumericalmethodfortwophaseflowwithanunstableinterface,J.Comput.Phys.39(1981)

179–200.

[16]J.Glimm,X.Li,Y.Liu,Z.Xu,N.Zhao,Conservativefronttrackingwithimprovedaccuracy,SIAMJ.Numer.Anal.39(2003)179–

200.

[17]R.S.Hankin,TheEulerequationsformultiphasecompressibleflowinconservationform:simulationofshock–bubbleinteractions,J.

Comput.Phys.172(2001)808–826.

[18]C.W.Hirt,B.D.Nichols,Volumeoffluid(VOF)methodforthedynamicsoffreeboundaries,J.Comput.Phys.39(1981)

201.

[19]C.W.Hirt,Volume-fractiontechniques:powerfultoolsforwindengineering,J.WindEng.Industr.Aerodyn.46–47(1993)327–

338.

[20]X.Y.Hu,B.C.Khoo,Aninterfaceinteractionmethodforcompressiblemultifluids,J.Comput.Phys.198(2004)35–.

[21]X.Y.Hu,B.C.Khoo,N.A.Adams,Asmoothedinterfaceinteractionmethodforcompressibleflow,in:The25thInternational

SymposiumonShockWaves,Bangalore,July17–21,India,2005.

[22]C.W.Huang,T.Y.Shih,Onthecomplexityofpoint-in-polygonalgorithms,Comput.Geosci.23(1997)109–118.

[23]A.Jameson,T.J.Baker,N.P.Weatherill,Calculationofinviscidtransonicflowoveracompleteaircraft,AIAAPaperNo.86-0103,

1986.

[24]G.S.Jiang,C.W.Shu,EfficientimplementationofweightedENOschemes,J.Comput.Phys.126(1996)202–228.

[25]S.Y.Kadioglu,M.Sussman,S.Osher,J.P.Wright,M.Kang,Asecondorderprimitivepreconditionerforsolvingallspeedmulti-phaseflows,J.Comput.Phys.209(2005)477–503.

[26]R.J.Leveque,NumericalMethodforConservationLaws,Birkha¨user,Basel,1991.

[27]T.G.Liu,B.C.Khoo,K.S.Yeo,Ghostfluidmethodforstrongshockimpactingonmaterialinterface,J.Comput.Phys.190(2003)

651–681.[28]R.Lo¨hner,Theefficientsimulationofstrongunsteadyflowsbyfiniteelementmethod,AIAAPaperNo.87-0555,1987.[29]H.Luo,J.D.Baum,R.Lo¨hner,Onthecomputationofmulti-materialflowsusingALEformulation,J.Comput.Phys.194(2004)

304–328.

[30]C.Mader,NumericalModelingofDetonation,UniversityofCaliforniaPress,LosAngeles,1979.

[31]M.Mandella,D.Bershader,Quantitativestudyofcompressiblevortices:generation,structureandinteractionwithairfoils,AIAA

PaperNo.87-0328,1987.

[32]G.H.Miller,P.Colella,Aconservativethree-dimensionalEulerianmethodforcoupledsolid–fluidshockcapturing,J.Comput.Phys.

183(2001)26–82.

[34]S.Osher,J.A.Sethain,Frontpropagatingwithcurvaturedependentspeed:algorithmbasedonHamilton–Jaccobiformulation,J.

Comput.Phys.79(1988)12–49.

[35]R.B.Pember,J.B.Bell,P.Colella,W.Y.Crutchfield,M.L.Welcome,AnadaptiveCartesiangridmethodforunsteadycompressible

flowinirregularregions,J.Comput.Phys.120(1995).

[36]J.Peraire,M.Vahdati,K.Morgan,O.C.Zienkiewicz,Adaptiveremeshingforcompressibleflowcomputations,J.Comput.Phys.72

(1987)449–466.

[37]A.R.Pishevar,AnALEmethodforcompressiblemultifluidflows:applicationtounderwaterexplosion,in:CFD2003,Vancouver,

May28–30,Cananda,2003.

[38]J.J.Quirk,Analternativetounstructuredgridforcomputinggasdynamicflowsaroundarbitrarilycomplextwo-dimensionalbodies,

Comput.Fluids23(1994)125–142.

[39]R.C.Ripley,D.R.Whitehouse,F.S.Lien,Effectofmeshtopologyonshockwaveloadingcomputations,in:CFD2003,Vancouver,

May28–30,Cananda,2003.

[40]J.A.Sethain,Evolution,implementation,andapplicationoflevelsetandfastmarchingmethodsforadvancingfronts,J.Comput.

Phys.169(2001)503–555.

[41]C.W.Shu,S.Osher,Efficientimplementationofessentiallynon-oscillatoryshock-capturingschemes,J.Comput.Phys.77(1988)439–

471.

[42]R.W.Smith,AUSM(ALE):ageometricallyconservativearbitraryLagrangian–Eulerianfluxsplittingscheme,J.Comput.Phys.150

(2001)268–286.

[43]M.Sussman,E.Fatemi,P.Smereka,S.Osher,Animprovedlevelsetmethodforincompressibletwo-phaseflows,Comput.Fluids27

(1998)663–680.

[44]E.F.Toro,RiemannSolversandNumericalMethodsforFluidDynamics:APracticalIntroduction,Springer,Berlin,NewYork,

1997.

[45]A.B.Wardlaw,H.U.Mair,Sphericalsolutionsofanunderwaterexplosionbubble,ShockVibr.5(1998)–102.

[46]T.Yabe,F.Xiao,T.Utsumi,Theconstrainedinterpolationprofilemethodformultiphaseanalysis,J.Comput.Phys.169(2001)556–

593.

[47]G.Yang,D.M.Causon,D.M.Ingram,R.Saunders,P.Batten,ACartesiancutcellmethodforcompressibleflows.PartA:static

bodyproblems,Aeronaut.J.101(1997)47–56.

578X.Y.Huetal./JournalofComputationalPhysics219(2006)553–578

[48]G.Yang,D.M.Causon,D.M.Ingram,R.Saunders,P.Batten,ACartesiancutcellmethodforcompressibleflows.PartB:moving

bodyproblems,Aeronaut.J.101(1997)57–65.

[49]D.L.Youngs,Time-dependentmulti-materialflowwithlargefluiddistortion,in:K.W.Morton,M.J.Baines(Eds.),Numerical

MethodsforFluidDynamics,AcademicPress,NewYork,1982,p.273.