An interface tracking method for hyperbolic systems of conservation laws

Applied NumericalMulhclnalnrsIO (1~)~12),147-472 North-Itolland

447

APNUM 351

An interface tracking method for hyperbolic systems of conservation laws Stephen F. Davis Ct~rh,R44, Naral Sllrfiu~' IVarfan' C('~ller. White Oak, MD 20903.5000, US.,I

.,lb~trac!

Da','is, S.F.. An interfacetrucking mcthtld for hyperbolic systemsof consev.;,tion la~'.s. Apl~licd Numerical Malhematics l0 (1992),147-472, This paper describes a method for tracking contact discontinuities fin(! ma! ~ial interfilees that zrlse in Ih¢ solution of hyperbolicsyslcmsof (2OllSt21~;itiolqI',lws.Numericalres,Jlts arc W :nlcdlil ~htlw ltl;11 Ihe fronts are rc.,,olvcdto v,'hhin a nlc~,ll int.w,'al and f;nloolh porlions of tile soluIiOh are eonlpuled m within Ihe accuracyof lhe ullderlyirlgnUllleric~ds,chenl¢,

!. Inlnlductinn It is well known that when second-order Godunov-type methods are applied to the equations of gas dynamics, they do an excellent job resolving shocks but do less well resolving contact discontinuities. This behavior can be explained by thc fact tMt shocks are "genuinely nonlinear" diseontinuilics while cont~bet discontinuities are "linearly degenerate". For the discussion which follows, it suffices to say that lhcse terms imph-' that the characteristics of the Euler equations converge at shucks but are parallel at contact discontinuities, Wh,,:n these equations are solved numcric~flly, the stecpcning of the solution caused by the convergence of the characteristics at a shock is balanced in part by Ihc dissipation needed to stabilize the method. As a rcsuit, the shock is ~Jppruxir,aated by a smooth but narrow transition. Oil the other hand, since there is no natural steepening o f the solution near a contact discontinuity, the numerical dissipation in tile snlutiou scheme causes the approximation of the conrad discontinuity to spread until either nonlinearities i~ the solution ut aonlincarities in tire numerical method stop the spreading. Recent ssork by Yang [35] examines a method which artificially steepens contacts. He shows thai this method reduces the amount theft numerical methods spread contact discontinuities but it is not able to eliminate this spreading altogether. This makes it difficult to sulve problems which contain interfaces separ;Iting different materials. If. for example, the two materials have different equations of st~tte, some rational procedure must be developed In determine an equation of slale for tire "mixture" of materials in the smeared interface. We are not aware of Corre.~pondcnce In: S.F. Davis,Code I/,44,NavalSurf~lceW~rfareCenlcr, While Oak. MI) 2(1903-51X~1.USA. 11168.9274/92/$05,(~1 ~,~ 1992- Elsevier Sciencel'ubli~,hcrsB,V, All right~;rer,en,¢d

4414

s,E Darts / Trackillg ]br hylu'~qlolic ,9'sit'nix

any general procedure for doing this. F'or this reason, this paper examines the alternative approaclt of tracking slmrp interfaces which separate distinct muterials, [:or tile present work, wc ilssunle that the material interfuce exists when the calculation begins. This means thut we do not attempt to find interfaces which are created during the course of tile computation. A futnre paper will address this rdoted but separate problem. Likewise, we do not attempt to track shocks. As stated above, we believe that the Godunov-type ulc~huds do an excellent job of resolving shocks, so there appears to be no reason to track them. Tllcrc may exist special problems for which it would be necessary to track it shock but, in gel~cral, it woukl be better to avoid doing so. We say this because shocks arc subject to many complex inlcnlctions with other shocks and witb material interfaces, and because it is far more difliicult to determine the velocity of a shock than that of it material interface. 13ccause many problems in cotnputntional physics arc concerned witlt tile dynamics of nmtcrial interfaces, many nttmcricnl methods have been devised to solve tltese problems. An excellent summary of this work is rite re':icw paper of Hyman [19], llcrc we outline the basic ideas behind Ihe !nnst commonly used methods and their shortcomings. Then we dc'~cribe the apgroach that we propose in this paper and explain why we believe it to be particularly robust. 'File most straightfurward way to treat it problem which contains an interface is to use a L~ grangian method [11-13]. This means that the equations of motion are written in such it way that the unknowns are the locations of discrete material volumes. Thus tbe location of the interface is part of the solution. Unfortunately, material cells often bccotne extremely distorted and numerical methods which attempt to approximate a solution on it distorted cell become less accurate. Eventually. cell bnundaries cross and the computu'.ion hills. Tilt:; problem can be alleviated somewhat by redistributing the grid at varinus times during the calculation. Unfor!u,, natcly, Ibis redislributiou can be Ihne-consumhtg and can introduce additional interpok|tion errors inlo tile computation. A second approach is to approximate the front by a set of marker points connected by an inlerpolant. The evolution of rite front is determined by solving a set of equntions for tile evolution of the mnrker points. Examples of this approach can be f~uad in [14.22,25]. There appear to be two difficulties with this approach. The first is that complex data struelures arc rcqoired to specify tile location and cunnectivily of tbe marker points. The second is that numerical errors can cause the iulerpolants to become entangled and the computation to fail. These problems arc more severe iu three dimensions than in two. We are unaware of the existence of it three-dimensional marker particle code. Another approach espnuscd, for example, in [3,5,18,23] tracks the vnlumc of one of II:e two muterials through fixed cell interfaces. This is done by introducing us a new vltriable the vtJl',Jlee fraelino of the tracked material, The pattern of volume fractions in a set of cells is used to rectlnstruct the interface, and n local Lagrangian calculation is cltrricd out to advance tile interface mid update ihe u~lume fractions. The present author has never tested this approach but i;ublished results show it to be quile robusL Critics conlphtin that rite resolution of these tnethods is limited to the cell size but could more be expected from any numerical computation? The ilewest approach Ihat we are aware of is Hartcn's subcell resolution method [16]. This utethnd reconstructs an interface by taking adwmlage of the fact that the numerical solution can be interpreted ;~s tile average value of a variable within a cell. In addition it ix assumed that a single interhlcc is contained in it particular cell and tfiat the solution is smooth in adjoining

S.F. Dat'is / Tt~rckingJbr IOT~erholicsy.~wms

449

cells. With this information it is possible to construct a discontinuous function which interpolates the solution wdues in adjoining cells and places die discontinuity so that the average of the fnnction matches the known cell averagc. A lucid Eulerian c~dculalion which accounts for tile movement of tile intcrfikce across cell boundaries is used to update cell averages. The method also contains a elevur scheme for locating cells which coatuill discontinuities, Hartcn carried out sample one-dimensional calculations in his paper uml the results were spectacular. tJnfortunately, an ~ttcmpt to extend the nrcthod to two space dimensions using splitting w~s not so successfnl [3(I], Thus far, a multi-dimensional version of the method has not been devised, The method that we propose in this paper hus some features which ~re similar to the volume tracking and the snbeell resolution method but contains two csscmkLI differences. First, we attempt to uncouple the problem of locating the interface from the problem of determining the conditions ~at the interface. To do this, we note thai the type of interface cousidere:_' here moves at !lie local velocity. Therefore, following Mulder et al. [21], we define an arbitrary smooth function g which satisfies the differential equation, ~g -~ ~ + n ' Vg = 0, (l.I) where u is the local velocity vector. If we know the initial wtlue of g at the interhlec, we need only locate thai value of g Jt any later time to locate tim interhtcc. This lerel set fiutctiou is eltoscn to be smooth so that equation (I.1) can be solvcd accurately. This idea is simihr to the volume tracking method except that the level set function has no physical significance; ;t is an arbitrary smootil fnncfinu whose wdues move at the local veloeit3,. The volume fractions, on the othcr hand, are not smooth. Thus it is necessary to perform a careful Lagrangian calculation to prevent numerical spreading and obtain :m accurate location fnr the interfitce. The present method can use the same Eulerian method to locale the interface as is used to solve the rest of tile problem. By contrast, the subeell resolution method uses the pattern of solution values to locate the intcrfacc. This presents difficulties in morc than onc space dimension. The second difference between the method proposed here and either tile volume tracking or subccll resolution method is in bow the flux through the boundaries of the cell containing the intcrface is computed and how the interface is moved from cell to cell. The othcr two methods construct a det;liled picture of the discontinuous solution in the cell cuntainiog the interface. This constructed solution is used to determine when the interface crosses the cell boundary and the flux through these Imur~dnl'ics. The method proposed here attempls to avoid such detailed snbgrid calculations. To do this, the solution values nn each side of the intcrE~ce are uncoupled from each oilier and die flux is determined in such a way that the lilaterial in a cell corresponds to that specified by tile level set filnclion. When the material in the cell clumges, the jump across the intcrf;~ce is added to the solution in the ceil. Local approximate I,~,ienlanu problems are solved to uncouple the solution ,'rofi;es on each side of file ia~crface and to dctermi,~e the jump across the interface. In the remainder of this papcr, we describc thi:; nlethod iu detail and apply it to problems in one space dimension to illustrate its beb;wior. In parlicuiar, we show that this method permits slrong shocks and expansinn waves to pass through the interface ~md that second-order accuracy is maintained away from tile interface. Ill -'ldditimt, we brielly describe an operator split version of this method [32] which can be applied to two- and Ihrce-dlmcnskmal problems

S.b~Dacis / TrackhlgJbr hyperbdic.~ystems

450

•rod present sample calculations. A detailed study of the multi-dimensit}nal rnetbod will be presented in a flJture paper.

2. Description of tile algorithm In this section we describe tile algorithm that was developed for this study. Tile basic idea is to make tile overall ntethod as simple and robust its possible. To tlds cod, a smooth level set function advectcd with the local fluid velocity is used to locate the interfoce. Smooth solutions can bc computed accurately and other fluid v~riables do not dircctly affect the interface location. This approach was first proposed by Muldcr et ill. [2l]. Our contribution to this work is to derive special proccdnres to determine the flux through the boundaries of cells containing interfaces and the solution witbin interface cells. Thcse procedures permit pressure waves to pass through cells bat prevent material from behind the interface from entering an adjoining cell until the interface itself is well within that cell. Once this happens, the solution within the cell is set to an accurate approximation of tile conditions behind the interface. The purpose of this section is to present a general framework for the construction of an intcrfilee trucking algorithm. Specific implementations of different parts of the algorithm such as numerical flux ibrmulac and approximate Riemann solutions arc dcscrib~:d in appendices. The re.'lson tbr Ibis is to separate the general framework from specific implcnaentations. Other imple:,lentations can be incolporated within the general framework to produce effective methods for solving interface problems. Later sections contain the results of numerical experiments wbich denlonslrute the performance of the algorithm. Rigorous annlysis is deferred to future papers or other authors.

2. I. Method fiJr trackhtg front Assume that the inilial location of the front is known. Constrncl a smuulh, function g(x) which is monotonic ill tile ncighburhoud of the front and has a zero at tile front location. A function which behaves like

g(x) = C ( x - x j , )

(2.1)

ill tile vicinity of tile fr~mt location x~, might I')u n good cboicc. A material interface moves with the local nmterial velocity. The mathematical st;ttemcnt that the values of the function g nlove with the local vck)eity is dmt g satisfies tl~c cquati~m + u .Vg = II,

(2.2)

J

where u is tile local velocity vector. If equation (2.2) is solved with initial conditions which behave like (2,1) in tile vicinity of the interfilcc, the locution of the interface will correspond to the locution of tbe zero of g at all times. A nlor¢ suphisticatcd lint cuncepluully similar appn)ach has been proposed by Usher and Selhiau [24] for tracking fi~mts that nmve with a speed that depemls un tbc geometry of the front. A version of the present mctlmd, where the front moves wilh the local velocity, was first propused by Muldcr ct al. 1211. In llmlliple space

S.F. Dat'is / Tracking]or hyperbolics.~!~wms

45 I

dimensions, the initial rabies of the fllnction g arc computed as the signed distance (positive on one side, negative on the other) from tile initial location of the front. In tiffs study, equation (2.2) is solved using a non-conservative version of the second-order method used for the field equations. This method is described in AppcrldiX A. Since u is continuous through u material interfilee and the initial conditions on g arc smooth, g will be at least continuous. This means that standard methods [31] can be used to obt~iin error estimates for g. Under these circum.qtanecs we would expect second-order accuracy for tile location of the interface. The cell containing the interface is located by the following simple linear search of It cells in one dimension. flag(n) ,- no for i - 1, 2 . . . . . n - 1, do if g ( i ) . g ( i + 1) ~
452

S.F. DarL~'/ 7)'~lckingfi,"h37Ji'rbolic!.~3'.~lem~

flux into tile cell 1o tile right of the tloumlaz3' ix computed us tile physical flux evahmtcd at tile right-hund cxtrupolatcd solution value nlh]us tile contribution of all right-moving wuves except tile inlerf;icc. Tile flux into the cell to the left of tile boumlury is conlputed as tile physical flux evaluated ut tile left-hand extnlpol+tted sulutitm vuhle plus the enntribution of u]] left-moving waves exeept tile interfnce. This procedure assures that 11o material fronl tile other side of the interfuee crosses the cell boundary, lustcad, the contribution to the flux on each side of the cell boumtury frorn the eharacterislic field eorrespondir~g to lht2 eolltael, and only that churaclerisLi.,: field, is deternlined from u same-side extrupc,lati,.n't. In the last subseetioiL we present u hettristie argument to support this procedure. If tile values of g in tile two cells adjoining tile boundary have the same sign, Ihen the flux throtlgh file boundary is conlpuled using tile same method us is used for nou-interh~ce cells, "Ik) illustrate this procedure, consider cxtrup,31ated solution vectors Uj. and Ur assumed to be m-vectors. Using either an ~pproximatc or an exact Ricnmnn solution, tile differencc Ur - UL can be split into contributions c(,rresponding to different waves. Th~lt is, U r - Uj. = ~ h

~U "~,

(2.3)

i

;shore each AU ~t~ is tile vect~lr jump in tile sc~luti,.m across tile wave which nluves with v¢!ocity a ~ ~. Assume that index k = / corrcsp~mds to the interface. Tiles1 if g has different signs in tile twcl cells adjoining a cell boundary, the numerical flux into Ihe cell to tile Jcfl o f tile cell I'~oulldary is Ft.(Ul., Ur) = f ( U i . ) +

~..

a~ ~-U('~1.

(2,4)

I.I, : ~Jtt ~ .: tl)

t.¢l

Simil;trly, tile numerical flux into the (:ell to the right of tile celi bt:,und;iry is Fi~(Ut.Ur)=f(Ur)-

.~.

,"~aU'~L

(2,5)

I/, : ,t q~ I~b} t.~'l

In practice tile solution is cumpulcd every~'herc usinLl !he same method ;rod then corrections arc ;idded to !lie cells ;ldjoining a cell botinduly c,lmtuiniug an interface. The ¢orreclinns ~lssoci:lted with tile numerical flux described ill [ll)] are presented in Appendix C. 23. Sohahm i'ah : hehhul hlteffuce 'The procedure described in the previous std)sectilnl assures Ihzl onl.v mulcri;al frotn (Ine side of afl inlerfucc crosses a cell boundary duriug a lin~c step, Waves cross cell I,ouudarics unhindered, hi all cases the maleriul in a cell ix idenlified by the sign of tile level set function g, If the sign of g shoukl clmng¢ within u lisle slep. it is nc¢cs~:uy !~ ch,ng¢ file material in tile cell. This is ucconlplishcd by adding the jump in file snlutinn /aU "~ corresponding Io lilt: interfuec to the solution v;due in Ih¢ c¢11. The sulutiolJ jmnp is comput .d by solving a Ricmann problern with one state correspunding to the solutinI~ in Ihe new il;t:riace con and the other st;de ;111exlrupolation 1o tile centcr of tile i~cw JlllCltac¢ cell from tile ,old interface: cell, ']'11,2ILet effect (if this prnccduve ix to oh:rage Ihc IIlateriul ill U cell ill discrete iumps. This avoids the crealinn of ntulli-fluid cells, Sucl~ cells rctluilc either models hlr mixtures of fluids or

S.F. Dm'L~ / Tnwkin~,'fi~r h~l~erhdil' ay.~tcm~

453

. .u-lift

j--1

j

j+l

j~;2

..... ~lt | UL " - - -

_.O'- ..... t!~t." _ ~).-~-r'- -ta R"

r j 1

I j

: j~l

: }~a

Fig. I. Iqam,ib ly atgumt:nl for the melht~t|. (a) Inilial condilh)rrsaml (h) condilion~ ~hcn lchrd NO1fllllClil!lltl;t~, challgcd sign. Da:,,hcd lines d,ant~tecxlrapt}];lled'¢:lhJ ¢ *i. speehll dat~l ~tructures Io sR~re multiple solution v~lues. Rdow we show that, fi~r (~ne-din~etlsitmal prohlems, it is still possible m obtain accurate estimates of 111¢c~mtael jump and Incat~(m. 2A. Ihm, it ww'ks

In this subseclion ss,e give a plausibility argument to show thai Ihe melhod described above exactly resolves a dixconlilltlity ~¢parilth)g l~.:o linear profiles. These arguments are suppt)rtcd by numerical experimcnl in a later seclion. To keep things simple, consider Ihe equation tqlt ;~lt ,'t7 + ~ = O, (2.6) The initial dala is a discontinuity at pninl j separating two linear ftlnclions with different slopes ;'.s shown iu Fig, I(a). Cell /' ecmtains the value n~ to the left of the discontinuity, u~ is the value on the right of the discontinuity, and r q atld it R are linear extrapolations Ii) the hotlndary hetwcen cells i arid ,/+ I from the left and right respectively. In this case, the procedure described above ~,,ill compute n n a~, the flux on the right ~,idc ~lf file cell boundary and u~. as the flux on the left ~,idc of the cell boundary, Thb, discon,etls tile solutkm profiles ~m each side of the discontinuity from each other, Since the mlutkm profiles are moving Ii) the right and the numerical methods used here :lle exact for linear profiles, Ihe exltapolatkm of the right.hand profile will propagate into cell j + I and the ex|rapolati~m of the left profile will propagate out of cell .L

4.~4

XF. P.aris / 7),leki~lgJbrh.ffle~bnlic.~y~l,'nt~

When the level set fullClion calcuhlti¢ln, carried ¢lUt ind¢l'~cnd¢lllly, btdieales Ill;it tile discontinuily has m,)ved Io point j 4- I, the suh=lion is as sll~wn ir~ Fig. ](h), Tilal hie;ins |hal the solutiem ill cell j + I is ,~. A linear extrapolation .,If the solul/o;l in cell j hl lh,.: cenler ¢if cell j 4- ] is tt~. There[ore fhe inlet[ace jt=nlp is u~ - u] ;rod lll~.' new solution J, cell j + I is u~. This is the resull we started with Iranslalcd hy one grid i~llerval, Note Ill;it v.,¢ have assumed lhaI the inferfitce nlovcd exa¢lly from point j IO poiltl j + ]. hl practice, lh¢ irll¢fface will not fall ~x;Iclly on a/:;rid poJlll, This Silould ~.'ause no Imflllelll fi~r (Hlc-dim~:nsional calcldations because .lie cmnputed juml) slill c~lm)ccts file fwo It.ear profiles. The e,'~a¢l ]Ue;lli~in can be dclernlined from the level sol fi~lelicm, and an ;ic.~'tlral¢ jlllllp ¢;III D~: delermi;led by c×Iral)olati~i[; from each side to this hlcafi~lll.

3, Tv, c~-dim~:nsio~lal e,~lension

hl lifts ,,cello;l, we extend the nlethod of tile preceding scclioll Io Ira)re ~hall o11c space dilncllSiml via spliitiilg [32]. 'ro apply splitting, wc ~Lssunlc |hal the proNenl domain has bccn c~wcred v,'ilh a logically r¢ctanDflar mesh which is paramclcrizcd by two ird~gers i an~.| j. Apply lit,,: algorilllm (ff lh¢ l;isl scctiol~ ;ll()rlg each j = COll$t. line to a set of ¢qtlalions ,¢:(lllStrUclcd by setting ~dl j deriv:~lives Io zero and dclelinlz source terms. I)¢llL)le this t~peralilm hy Ihe Synlb(lJ L,. I.ikewise, dcnolc hy Lp tile al~plicaliotl of file algorithm ahlng i ~: CO~lSl. lincs to equations co;lslruelcd by dclclhlg source terms and terrns involving i derivatives. Lcl L~ ~Jenote Ihe solulk~ll ,a~ Ih¢ equati~ns constructed hy deleting all terms involving b,,llh i and j derivatives. Two lin|¢ . :¢ps o[ Ih¢ .';,plitting scheme are dclloled sy1~lh~flicallyby w"' : o: I.,LjL , L,I.~L,w".

(3.1)

This nlc;ms lhal lh¢ i ~ipelafi~)n is aplflicd h~ the hlifial data, file j opcrali~ii~ i~ applied hl tilt rcsull, file ~,' ~'~perali~n is applied fo lifts re',,ull and tlivll lhe OF,clali~,llS me a|;.plied ill turn hl this resuh ill reverse order. Slr;ing [32] Ii;l~; shown lh;ll if ¢;ich of tile COnlpOl;enl OpCl;lI(,~l~; .'Ire second.order ;ic¢iir;ii¢, then t llc ~plilling SCllCnle (3.1) will he sceond-=irder ;Iceur;ll¢. The Iracking scheme w(;uld be idc,lical 1,1 |hal used ill o11¢ ~,p;ICC dilncnsion if all eli' the: linearly degenerate W;lV¢~, were contact disCOlllinllitJ¢~. Ullf~rltJ1~;~l¢ly. ill tWO ,lr more ,~p;~c¢ djnlCnSjOllS, lll¢: .~;=~dyllalnics cqualion'., pos,,c,~s ;u~ol]le: Jhlcally dCgCl|Cr;de w;iv¢, lhc sbc;ir wave. Shear waves move with the sanl¢ ~,'c]o,'Jfy ;is cllnl;icl diseonlinuilies hut lhc~/ are ¢lricllt;fliem-dcpcMenl ;llld thclcfllle ;Ire more difficult h~ resolve numericall~. 'l'hcse dilliclJl. tics are e'q)lained heh~w. In this paper we avoid Ihc pl,~hlen~s as~,~Ci;lled Wilh she;=r waves hy COllSlrIt¢lir~ ;i nIcthod wllich only tl;lgk~; collt;l¢l dJ~,colllillUilJch, This tile;IllS Ill;**'. file ¢llle-di. nlensi(~,;d ,llx.'ralors ;ice nlodified S,l |hal ,,,h¢;=rv,';ivcs are c;iplulcd and aIMwcd f=l spread. We ;ire CUlleBtly v,'oikin.~ tel d~.'veh~p 1111oriCllf;JtiOll-dep~...ndelll nlcfhod which ll;leks ~,hear w;ives, 3.1. Shear

wa!v.~ , n d c~mt,'.ct*'

3"o clarify i~ue~, we define :s ~hc;ir wave ;1~ ;l dkconliI~hy of the 1;UlgClzlial velocily COfllpOItL'rlt%, N()II|I~II vcloCl|y cofllp~lneflts ;1ic Ollllhlm~u~, This 111¢a11~thai lh¢ IIliCnlllliOn ()l the wave i',; dcI¢llnined hy the jump ill lhe velocity VCCt~II, II we km~w c~in,liti~in'., al twt~ l',,ainls,

,,ve :l,..'ed fo tl¢l~2ffllill¢ al'l orientation in order Io fit a shear wave between the point~L 'Ntis ',s not c~lsy |o do 11¢¢31:,'L'l[Ig SO[Ilti(HI jumps ;LCF()SS:lOrllill¢~lr pressure wa~¢s also depend upon ()rie:|lation. II" fh¢ (}rienlat[,,)ll CllOsen c,orrespollds Io a parlic~dar w~,ve, the re~olu|ion o[" that w~Jve is improved. This was shown in [9]. It" 1[1¢ wrong orienlalion is chosen, tile mlulhm to the [~,ienlarlll problem contains nlUllir~le ~.vi:ves. T[|[~; n|CilllS Ihat a soluligm wl'&h conlilills only a shear wave if file oriclltafion is corr¢',:l, 13)llklill~, [wo pressure waves alid a weaker shear wave when tile oriclllaliOll is :loL corrccL Lil..:ewis¢, a pressure wave wtluld be interpreted as a shear wave ~md two dilfclcn! f, rcssule W;lyeS [[" Ihe or[¢nla|i():l were not eorred. This dependence on oriental:on cavses coordinate split versions of mclhods such as II]¢ random ¢hoi¢¢ mclhod, which cxaclly resolve individual waves, to prodLi¢¢ kffg¢ errors near waves Iha¢ arc :lot aligned will: a C.~.)ordillat¢ direclion, This problem wa:-; analyzed in [2~I]. This aualysis explains wily il is difficull Io conslrucl a splil nlethod ~,,~lliell mleks these waves. By conlrasl, a ConlilE[ disconlinuily is ~l jump in a se~dar v~Lriab]~. |n gas dynamics il is lh~.~ d¢i~sily, This means Ihat lhc jump in Ihc solulion is independent (,f Ih¢ oricnlalion ~f Ih¢ wave and il is possible Io fit a eonla¢( hetween any two poims. Because af lifts, we can Iraek Ihes¢ w~lves with a coordinale split numerical m¢lhod. A nl;llhematieal cxpianafioii ,if dli~ behavi~}r can bc i[luslral(:'d by sludying Ih¢ El:h:r ,.qualions in two space dimcnskms, a.

af(.)

a~(:,)

+ "~-~-'x + ~

~ "'

(3.2)

~'here

. '.~ l~,, ira, p,,, E] I ,

(3.3)

f(.)" It,.. ~,.= +~,, ~,.,', (r +~0.]', ~(") '" l~"'. p"". p"~ + ~', (I; +,),,]~.

(3,~) (3.s)

For simplicity, we consider lh¢ perf¢c| gas ¢(ItlaliOll o f ~,(al¢,

~, ,~ (~ - I)[ Z; - ~,,(.: + ,,')].

(3.6)

The cigenvalucs and ¢(irrcsponding cigenved~)rs of Ih¢ .lacol;iall irl;llliX A ,* [;tJ/atl] :Ire

~l~u+a, A[~:,-a,

r~,[I,u-la,:,,ll.l-ua] ~, r~-~ [1, I ¢ - , L v, It-:m] ~,

~', ~., ~: °~ 10, o. ~. ,,1~, Ily symmetry, lh¢ eigenvalues and eoeresp,mding ei.vew,'¢clors of II ~ I;~/~/au] arc

~=:,-r/,

A]~v-l.a,

r~[1, u~v+,,ll+;,u} t, r~[1, u,:,-a, lt-;,a] ~.

~:=,,.

,~, = I~L i.,,,,] r .

S I" D,.i~ / 7~. ki~lgJbr h~/wffr~lic .~l~t~vn~

451~ In Ihc III)|IVC

11 = ( E

+lOl:',

(3.~)

q: ~ I¢" + v ~,

(3.111)

a: = ( ' / - 1)(11 - ':q?).

(3.11)

The cigcnv¢clors r~ and r,~ corrcsporld Io eOlll;icIs and shear w~wcs respectively. The other I~;o ¢igenvct.l,t~l~ o.lrrespolld to pressure waves, Nolo that Ih¢ ¢igenveclor r 3 is the only one :l,;it is ct.umon t~ both illiltli¢cs _.1 and B. To ~ce Ihc significance tff this, assume tllat the sllace derivatives t f f . arc i'lrop(irti(:qlal t~l this ¢igcnvcctor. That is, --r~,

- .....

(3.12)

~ll~crc. i,i a sc[d:lr funclioll, If these expressions are substituted into (3.21. fl~e rcsull is a single, uncl~lll}l(2(l SC;II;IFCqtl[lliO~l li~r t~. If v,'¢ ;J~;,SUllleIhat Ihc Sl'~;Ice dcriv;~tivcs of is are i')l'optwtional to r~, the equations will not IlllC~,ltll'~l(: b¢caiP,¢ r~ is Ilo! all ci!.cflvccttlr of B. II is c:lsy to show Ihal Ibis vector will project ,tllllo Iht: three eigcnvcctors r I. r~, and r4'. This is whal wc tncan when wc say IllaL lJ shear w~ve ill ()lie dirccliOll is inlcfprcled as three w;ives ill c',.'cl'y (~ther direction.

3.2, Algorithm cha.ge~ As p~linlcd (~ut ill Ih¢ hlsl so(lion, the IIonl h)¢~tioll mclhod ,~lf Muldcr el ;d. [21] carries ~wer to fflullil'~l¢ Sl');~,:e dinleflsiii!lS with csscnlhdly II() change, hlilial colldiliol)s tilt tile level set ftlll¢tioll .~, are d¢l'i;~cd by ;~,signing. Ill,,2 ~'13'o ,.'onlour Io Ill,,: initi;d Ioc;]lion of Ihc il',lerh~¢¢ and ddining the vahte (ff g to bc Ihc sigiIcd distance from lift',; zero C(lllt,t)tll. Since we allow .',hear i;lyer., It) ~,pl¢;ltl. Ih¢ fhJid ,,.,chwily u is i:a~lllilltl()tl% Ihrough tile interf:~ce. "l'llis means thai Ihe eqll;lli()ll ~,atisficd by g. i}g ;~t "~ u ' Ve ~' t).

(3.13)

h;ts conlinu(m., coell'icicnts and c;m bt! ~..(ll~cd by ~plilling, The ,mh,,:r two p;Jrl'., of li)c algolilhnl ilp,,tll~.,¢ ~pct;ial (fc;l|lll¢llt J()r the ll;lll (if tile ,'s(I]tlli()n jump ~,thieh COll¢~pOll!,J~, |o t" illl¢lf;l¢t; !n l~l)l)l~. |l~,ll) Olle ~0;ICC dhllcllsi,.lil it is illlp,.~rlanl h) ;,l~l'ly these pfoccdul¢,, 1(i Ill...tHll;icl linty, not lilt sl~c;~r wave. For g:ls dynamics this Ille;ln.';, Ill;l| special pl(~cedu~cs arc ;ipplk'd to ll'm vcq[ol, aU"

::

l ~,," -- at," l~:] ~,.

(3.14)

WJlcrc tile Iclm'~ in (3,14) ;,~c evaluated dcpelld'.i Oll file I~Hlieular :~ppmxin~ate I~.icmann

'..llution, 1:(~I the rn¢il)od described irl Appc;Idix B, Ihe "l'~ar" lernl'~ ale cwduifled ~t the olntzt¢l, ~, ill /'~ i", I;tken It) be tile ;Hhhnl¢lJ¢ ~vcr:,gc of the conditions ~]cros,, tile .,,hear lay(r, ~ p " = I., :nnd .'~p" is Ihe dew.ily dilievcncc across lh¢ ¢ont;]cl. Fer the n~clhod of R~e 127], all " b a r " terms ~ e " R o e ;~vcragcs" and ;dl ~ lernl~ are differed)cos between left and right solution '-,l;]t¢'~,

s.F. Din.iv / TracMng[or hYlll,rbolirsystems

457

4. Numerical results Ill this seetinn wc present tile results of u nuulbcr of numerical experiments to illustrate the performance tlf the method. In particular, we show that the'method is second-order accurate in smooth regions of tile solution, first-order accurate near a discontinuous interface, and that a smooth region treated as an interface does not lose accuracy. Also, we show that the method allows pressure waves (i.e. shocks and expansions) to pass through an interface without loss of accuracy and the method handles interfaces between vastly different materials easily. Finally, we show that this method can be used with splitting to solving multi-dimensional problems. 4. I. lixample I

For this ex~tmple and the one following, consider tile Euler equations of gas dynamics with an ideal gas equation nf state. That is, Ou --+ ill

a/(u) ~.r

=0

(4.1)

with "

= [o, p", el 1,

f(n) = [p,,, I,,," ÷ t,. (e +/,).] ~,

(4.2)

and p = (~, - l ) ( e -

½pnZ).

(4.3)

The purpose of this example is to study tile accuracy of the rnethod. To this end, we conslruct an initial density Po and an initial level set function go on a periodic domain as follnws. For tile initial density construct a Hermite cubic that interpolates p = 2.0 and p = 0.4 with zero slope on each end and displace it so that the discontinuity is in the center of the dom~dn. The level set function get is taken to be the union of two quadratics. The first, defined on the left h;df of Ibc tlonlaJn, is concave and interpolates ze¢o at the left boundar~ and the center of the domuin. The second, dcfincd on the right half of the domain, is convex and intcrpolalcs zcm at the center of the domain and on the right boundary. The solid lines on Fig. 2 show these initial conditions. The rcnlaJningpar:un¢'ters a r c taken to he consu.nll. Tile results shown use u = 1 1),

p = 2.0,

"r = 1.4,

Ar/Ax ~ 0.~.

(4.4)

Wilh tbcse choices the fl,aw is subsonic and the contact discontinuity moves through a cell in four time steps. In the e~act solution, these parameters remain constant and the original denfity profile translates with velocity u. The density sulution could have been studied by solving a scalar linear advcctinn cqualinn but the entire system was solved in order It) determine if coors in the method v,,ouhl produce spurious pressure waves. Fortunately, they did uot.

S.I". Dot'il / 7?atking fin" ll>pl'rhoh'c,~yizelm

458

0()15 I

,

:.

,

, ,

~

,

~- ,

i

I

() () 1 o [

i

MI

i !

i) [1

11.$~

(14

Og

0 II

I O

() 0

(I ,2

×

0 4

.....

0 (i

0 tl

I/)

X

Fig, 2. Solution tlf Example I with N= 50 after 196 lime steps and A = IL25.

Fig. 3. Local error in dcn,.ily in Ex;imph: 1 ~iller liJ(i lime steps ¢llmpulcd ',lilh N = 50 and J = 0.25.

Table 1 shows the eiror in Ihe density after it had returned to its original location, This means that the number of time steps in each case was equal to four times the number of space intervals N. The cell containing Ihe contact discontinuity was included in the error calculations to see how well the discontinuity was prese~'ed. The order of accuracy shown in Table I was determined cxperiolcnlally by Richards~m extrapolation. This shows that the computed solution is second-order accurate in the L= norm and slightly better than first-order accurate in the L~ norm. The maximun~ error occurred in tile cell containing the contact discontinuity. A lypical result for the density and level set function is shown in Fig. 2. The corresponding local error in the density is shown in Fig. 3. Note from the figures thai the level set funclk:.n crosses zero a second lime in a smooth region of the density profile. This =neans that the algorithm treats a ecll in this region as if it contained a discontinuity. The figures show thai this causes a slight loss of accuracy at this point but the convergence appears to be second-order ;is the ulesh is refined. No loss of

Table I I-~rlor i n den,,ily fk~r F.xample I after d i , , c ~ m i n u i t y r¢luin,, Io o l i r i n a l I o c a l k l n

] f11Cl'v~ll~,

L I t'fr~tr

Older

Ma~ ¢lior

Order

]3 25 50 {(~l

2.49555( - 2t 5.75036( - 3) 1.32749( - 3) 3.147fd5(-4)

2,24 2,11 2u7

9.22329( - 2) 3,94695( - 2) 1,4~33g( - 2) 5.2978l(-3)

1.29 1,41 t.4~

200

7A8973( - 5 )

2.07

1.8.! 24tK - 3)

1.53

400

t.743~1( - 5)

2.|(I

6,l 1(~97(-4)

1.58

459

.c F. Dat'is / T,':J('kingfor 107n'dm!ic,Lvst,'ms

~I

r.

,rl,

,r I ....

r,

.......

b

1) |

0

, ,r¢

20

,m

(J(]

II0

I OO

O

~11

>:

40

60

;¢'#1/ no

1oo

x

Fig, 4. Dcnsit:,,and le',cl set funclinn solution In the LIX nloblcrn ;lllci' 85 time steps wilh a C|)nl'alll number of I),8.

Fig. 5. Dcnsilyand t¢;'cl sut function,,,olufionto Exampie 2 alter 2(10 lime slopS. The shock has rcllected ffl~nl Ihc wall lind passed through the cnnl;ic| dlsconlinuity.

smoollmess is visible in the computed solution, even on coart~e grids. An earlier version of the algorithm introduced an artificial O(h) discontinuity in this region. 4,2, F~xample 2 The next test problem is ~l Riemann problem with initial data [p.,]

l;;:i]

I-(I.445]

=/°'6'~s/'

[p.]

[o.5]

[",,'~/=/°'° /. ,,-,,.,

,4,,

Tile level set function g was taken to be linear with value zero on the cell boundary containing the initiul discontinuity. Absorbing boundary conditions are imposed at the left boundary und the reflecting wall boundary conditions ate imposed al the right boundary, ltarten [15] attributes this problent to L~lx. The most prominent feature of its solution is a I~rge contact discontinuity. Figure 4 shows the solution after 85 lime sleps computed with lt30 points at u Courant ntunber of 11.8. The solid line is tile ex~'~etsolution and the symbols are computed v~lhles, Note that the break-up of tile initkd discontinuity into waves has distorted the level set function and omdc it only cc.nliltuoos. A careful examination of the solulion shows that tire computed contael lags the exact solution by ;JpproximMely one tonsil interval We suspect that this is duc to the fact that it takes a number of time steps for the computed contact to teach its final speed. The cx~tet contact reached this speed instantaneously. This error would not occur flit an isolated ¢ont;lcl. Despite Ibis starting el,or, all waves ale rcsol',ed ;tecvrately :lnd move at the proper speed,

4611

S, F, t )¢u'i~ / Tr¢lckingfor h.~7~l'rholic.~ystems

(a)

,

:,

,rrrl,,r,'',

or.,

't,+. =,,

(b)

'1 ....

1 ....

1 ....

i ....

i

5

?

,1

:1

,y, f,a,

..... ] .... ] .... [ .... I.... 11111 ~O 4(I 60 no I(J(J x x Fig. ft. (a) Den'dly;intl le'¢elset [unelilltl;intl (b) prct,sure soltuilm to Example 2 after 2411lime t,leps. An exp~tnsion ~'~l~,'t~is p;ISSillglhr(Ittgh111¢cmllilcldiscofllilluily. '

i

en

,I(1

¢"/]

i ,,

J I,,

Gn

, ,

ft()

Figure 5 shows the computed density nfter 200 steps. By this time the shock has reflected from the wall and passed through tile contact discontinuity. The interaction of the shock wiflt the conlaci has clcaled a reflected cxpansiou which is moving to the right towards the wall. This is the expected behavior but an exact solution is not available. Fimdly, Fig, 6 shows the computed density and pressure aDer 240 xteps. The expansion wave. formed when tile shock passed through the contact, has reflected from the wall and is passing through tile contact. Once again, an exact solution is not available but the filrnl Of the solutiou is as expected. Note IIle change in tile slope of the pressure ~lt the conl~lcl. 4.3. Example 3

Tile uext example, an underwater explosion of a gas sphere, actual]y requires some type (ff fiont mlcking fl~r its solution. Varialitms of this problem have been solved, using different metlmds, by 13allhaus and Holt [2], Li and llolt [20], and Charfier and Tessieras [4]. The equations of nlolinn are difli:renl in the gas bubbles ;Ind the surrounding water. Both sets of equations have the f,'~rm

~l, --

ill

;q(,) + . . . . .

ilr

,,(r.

,).

(4.6)

where, in the gas. n, f. and the equation of state are given by equations (4.2), (4.3), and 21¢ ,' = 7([,.

,,.

e +/q'.

(4.7)

S,F. DarL~ / 7)'ackingJbr hlTwrbo&.,~ystcnts

461

In the water° " = to, tm] 'r,

f(u)

: [+,,, ;,,,- +;,]

T

,

(4.8)

2u

,v = ~. [p, [,t,] "r, and p =p~" - ].

(4.9)

Equation (4.9) is a sealed version of the Tait equation of state [6,8]. The constant ?~,¢is taken to bc 7.IL Note that we have more variables in the gas than in the water. This causes no problem

for our procedure because the solution to the Rienmnn problem that specifics conditions at the intcrfi~ee is formulated in terms of the density, velocity, and pressure not the total energy. The particular problem considered has initial conditions:

I,,=, [o] [0] I Q | = i [)'0 /'

u+ = 0.0 .

pcJ L2.75]

I,~

(4.10.

0.0

The initial interface is located at r = 0.3. It is important that tbe water and gas be separated. The Tait equation of state implies that the speed of sound in water remains real even for slighlly negative wdues of the pressure. Negative pres.~ures in the gas would caase the calculation to fail because the speed of so~,nd would become imaginary. Figure 7 shows the solution to this problem after 85 thne steps computed with 100 paints at a Courant number of 0.75. Note that the interface is perfectly resolved in the density plot. The pressure aud velocity change slope at the interfitee and the solution contains no spurious oscillations. Compare these results with those in [4]. Figure 8 shows the density for I;::'. same problem computed on a 5(1 x 50 grid using a time split [32,34] version of the method. If the ordering of the split steps is reversed, the slight asynunclry in the contours shown in Fig. 8(a) moves from the x-axis to the y-axis, Tiffs indicates thai the slight lack of symmetry in these results is due to splitting error. Figure 8(b) is a three-dimensional plot of the density superimposed on t;le computing grid. This picture clearly shows that the interface is resolved to within a single cell in this two-dimensional computation. 4.4. lS'anlple 4 The final test problent is lhe computation of the unstable motion which occurs when a shock passes through an interface between fluids with different densities. Riehlmyer [26] constrncted a theory for tile initial growth of tile instability, and Meshk(w [1] studied this flow experimentally in a shock tulle. This is wily this process is often referred to as the li'.ichlmyer-Meshkov instability. Yuungs 136] has carried out computations to simulate the instability using a different numerical method from the one described here, Our calculations are performed on a non-dimensional version of a case computed 5)' Youngs. In particular, a long channel, 5 units wide and dosed at tree end is filled with a light gas (call it helium) from the closed end to a distance of 16.9 units. The rest of the channel is filled with a

S,I ,+, D,++'~v

/

7~.ckiPIl+

J?n'lijT~e+ hnli 08

; ~.y,+t,'m.+

(~)' l~ :

,;

I

"1

08

, ~?'+++~+

..... ~.+.+:-'-'+-Z-'..,=._~

04

g%,,

08

0

, ,e'" O(

o,,o~' 0.~

(1 11

112

0.'I

1}(~

(].ll

t) II

IO

{12

[1 ~

1] lJ

11 ~

1 11

x

X

:~ =~
O (I

11:2

Cl4

{I(i f) ~L ]O x Fig, 7. (a) Dcnsib' and tcs'ctset l'unctim],(h) vchlcity, and (c) pl¢ssuI¢ s
more dense .~as(air) :rod [he irllcll'accbClWCCll the ~ascs is pcrtt,rbcd hy ;Ill~Irlli}titlt ~'given by lhc equation ¢=tfl~ cos(l).4"~y), O.O <), ~5.0.

(4.11)

The amplitude ~F tile inili~ll disltlrlmncc, ~z.. is v;lried. At file initi~d time an M = 1.3 ~h~lck is prop:ignited Ihn;ugh the ;fir Iow;irds tl;q." h]|cl'EJcc ;rod ch:scd end t~{" Ihc channcl. Thi~ gcomctry is skclchcd in Fig. 9 ;rod |he in~,i;~l co'ndi|ion~ ;Jrc specified in T;lhle 3, The slmck passes Ihmutdl Ih¢ interface, [dlccls at [he wall, ;rod Im~Scs

463

S.I': Dat'i~ / Tr.cking for Itylwrholk"s)'.~tl'm~' 50 40

~,/t/~,~

30

Fig. 8, l)¢n.,,ily rc',uhs ft~r Example 3: (a) contoured al]d (b) superimposed on computing ~rid. 1'he time split ;dg(Irilhlll','.aS tl~¢d~l[h]axial ~,yllllllClr)' ;i[I()LI111112 y-;IXJS',',';IS;l~xun|ctl,

16.9

Fig, tl. Gel)mollyelf E~;:~nlpie 4,

through the interface from the olher dircclion. Computations of these early stages are shown in Youngs' [36] paper. Every 1i111¢ a .shock inleracls will] Ih¢ interface, I',','o shocks are •)rnled. Oi!e is Iran~,nlilled

through the interface and passed out of Ihc domain. The other rcllccts from the it|lcl[;J~e alld moves towards Ihe wall. This shock is weakcl than lhe one which originally reached the interface. It reflects al the wall and Ihe process is repeated. II is mlt surprising that no exact solulions cxJsl fi~r this c,;mplex process. Richlmycr [26] c(mlputes the initial growth rate ~tf a sinusoid;d interfi~ce interacting with a single shock. N',) resulls exist for mulliple .qmcks or hh';2 time behavior. The c~,nputalions ~lf Richlm'~cr [26] Table 2 t~lili:d eondifiem', fi~r Example 4 V;Jliahl¢

p u l' p y x..,h

I.~II 1.51h 0.523 n,(J 1,8fl5 1,4

4LL~

('¢nler ],0 (LO o,n I.(I 1,4

4t~475

l,~i,~ht fl.14(I.~ 0,0 n.~] I,q l ,(,7

6~1,375

(a)

.

L •15

' '

50

55

.../

(a)

;]o

k

,+5

60

.

'I

. 50

55

~,0

(b)

~5

5O

55

Ib)

co

.~5

+,5

55

00 (¢)

Ic} 20

~0 ;',

50

++++

+0

imd ((') ~¢h'il~, '~¢khlls I11[ I('~;Itll]11¢ 4 ;+I T ~ ~() oil it Imc l!Iid. Inili:d pL'I(Lllllitli:lll IlII ~+I I.

-:'J

~i0

~J5

C0

I:iit. II. (~d I~'~.~:l set ct+r1!oUl'~, (hi dcllSily l.'4)IIl{)Llt's, told 1~.')',~.'h+dl'.' +'~.'i'ItII~ JilF l~,:+rtll+Ik" 4. ;it T ~ 511 till :I

i'II;llS~ rlid. tniti:tl ll~.'IttJlhillitlll +;ll " ~J+.

sh~+w IIud lh¢ mili;tlint¢ll~¢¢ i:lo~slhr~Jt¢depends on ;illof lhc p:iramclcrs of lhc pFohlcnl. "l'hc'~c im.ludc tl~c inili;d ;mlplilud,: ~md fr~.'qucncy (fl+ lhe inlctf;tc¢ pcHud~:Jtitm, Ihc dt'l,sily Htit) ~icrt+,,,, the inlcrla¢~.', mid lh¢ vclocily imluccd by lhe ~ht~ck. The sh~:,ck i~ m~ inlinit¢ly thin di',omtinuity ;rod tl~¢ illlCl;J~.-litfll I;l~,,:~ pl,~¢~: ill'~I;llll;Lnct:,ll~,ly. F(~r I]I~ ittl111¢liCal l'¢SuIts lhill flflhlw+ lhc ",ht+tk i', ',pt~.';td u','~r ;J numb~I t~l" i'.rid i:ll~.*l'v[II~ :Jnd lh¢ l~llio c)I" lh¢ inili:il p¢llulb:lli(Hl ilI11plillILI¢ Ill llI~.' i~dd il11¢l",';d f~¢ILl:,:,~, It:, h12 ;Ir~ in:,pOll;llll r, IIIIII:,ICICF. h is II()I ¢1¢;II" ll()x',' 1h¢ 11LJ111¢Iic:~l IC'~ulls lchH¢ I(~ I11¢ phy~,i¢;d pincer,',. 'I'hc~I~il¢. ',.vc 111¢~.¢III lhcs¢ IC~Ull'., I(I ~huw lh;ll 111¢ llll111~LliC;II 111~.'lhI:,tI i'~ :Ibh.' it:, I¢~,~11X'¢ ~'Ol1:,l:,]i¢;llcd ilM¢lli~CCS ~.)11 l'¢I;Ifiv¢ly C();IF:.,¢ ~'t+Inplilil1~ +'~tids. A ¢~ll~.'h]l sludy needs Io be dt)l;¢ ltl dctt+l'nlill¢ ',A hat phy,,ic~ is ;J~.'Itl:+lly b,.+illlz ,,mmlalcd. Scc +d,,o lh~.. I¢I:~I¢d discus,,i
Fi~iil'C|ll xllt)%~,'S the lC~ults+:iIIcI',~()l:,t)tl.dill1¢l:,SiOlIIl] IiIHL'llnils+(~,i:i £1)Illpllt;Iti()I1 111:1d¢ wilh 411 I~'1idp~inls in tile y-dilcctim|(~y ,- II,125):rod XI) p~inls within the oriFimll helium let!ion (..~.1 + (L21125L 'lh¢ inili;~l ;Unl')litudc <,f 1h¢ p¢ilull);~li{)ll w:t',, ,,¢I cgual so ~ JJ ili1d ih¢ ('(iiii:i111 lllllllbl~l wit., ~,~I h+ 11..~. It Ii){+k 1734 lin1¢ ~,h.'p'~ I(:, le;i(.'h ;i llt)Ivdil:,ICll~.it:,llitl Ii111¢ (~IL .~1). The i'iI~Ul¢ include,, k'(Hl|tltll+ ()1' the I¢','CI ',¢I iuncli(~ll, den',if.',' ciiIlt+!lllS+ ;rod x'¢lc, city '¢L'~'I(11~,. l:it~'Ul~.' II ~.11<)w~ lll~.~ re'gill",, ;I[Ict' ~(I llOlt-dil:,1~ll~i()ll;II Iiii:,¢ IIllit% (~f:l k'l)tllp!J|;lliOrl on II 41) )<. ~II i.~lid ( ~ .~ - {L422.~, ~ y - lI.2.~) if:, lhc ,:,r if~ir~;~l helium ten,iota AI.~;dn. we ,,¢I ~ , - ~.~ ;rod |hl~ ('(1111";Itll HIIl:,lhCl ttl (),,~, It ltl~+k ~0() limi." ~h'p ~, It) lc;+ch ;I mm.diln~_'n'~iolml lime ~II .~0, The s;Jm¢ Cl)lIdili:lll,, 'w¢I¢ ll~L'd I~I ohl:iill I]I¢ l'¢~uIls sh~}~.n ill FiI~. 12 ;~s in Fil;. II ~'xccpl lh;d w.c ',¢I ,'t~ ,: II,~,l. Thi',, t11~.';tns th;ll 111¢ inili;~l ;llllplilud¢ ol Ili¢ p¢lllJrb;fliOll 'wa',, 111¢ ~,;i111¢ :+,, I~i' I if!.. ]I). 'lh¢ inili;]l ;m~plitud¢ li,' l-il~. I I w;~ ~,pl¢;Jd ovc~ the '~amc numbcl ol ~rid l)(~inl~, ;~,~ kll l-i!!.. I0. Ohviously. lh¢ Fit. 1 ] rc'~ulls ~]f¢ (:I~)'~¢i' t~) lhe Fly., IIi Icmdi'~ llmi| ;]r¢ lh¢ Fir. 12 l¢'.,ull'~. -l'hi,, me:m,, th;ll the l:flio ,%1~ is ;ill irllp~ll;llll pal'amCf¢l lof III¢'~¢ cl]li'ulali(ms.

XF. l),ris / ?~'aJ,htg fnr ll)71e~hnhc.~y.s+,v.s

465

(a)

(a) ,

30

'°L

.... 45

,

5O

!' ' ,!

i~

i

+ 5.5

:

10

:

45

CO

,,, !

5O

.... :"<;',

lj

;/ "1 , I : " 1 i I 55

60 (b)

(b)

45

bO

b5

L © 4b

GO

50

55

.~Lo 040

40 30° 10

i

,

45 t,O ++5 t,O Fig, 12. (a) I£',¢1 ,~cl ¢OIlIOUIs,(h) dct~ily ¢OlltOtll~.+ and (c) x,.'locity ;e,:lUlS fief I~',,ample 4 .ut T " 5ll cm :~ ('(131"*.Cgild, hfillal pcltulh;Hiol) +1,- (I,.~1J, ~ AJ I'

~:0

It)

Ic)

"':....

'~5 50 55 60 Fi,,z. 13. (a) Ix.'vcl gel Ct~llttltlr~,, (b) dcn',ily ccmtout~,, :ind (+) 'cchicily"¢+~ClOl'~+oi l+'Aaml+l,:4 ill I " ~0 O11:| fineglid.Initi;,I llCl'lUlbillioll +t. ~ (I.SAJ t ,

Figure ]3 ~ihcws Ih¢ r¢',,ulls, aflcr 50 (im¢ ullJts, ()~ ;! ¢ompulaliOrl usil'lg Ihc ~qtn1¢ grid 311d Cotlr;inl numhcr as wa~ used Io ohl;fin Ih¢ Fig. H) t¢~,ulls, lierc IIz¢ iliiiJ;ll pclluzbalhm ampiilude bl (I.5~.x, Nol¢ Ih;it Ih¢ rcsulls in I:ig. 12 :Ire nol :t,i clo~,¢ to Ih¢ I¢~,u11'~ill Fi~',, I,t as IIIo~,¢ in I:i#. I ] ; , ¢ Io those ilz Fig. I0. This indic;ltc~ timt 11~¢ grid m,ed was Ic~:~ o~;tr~,¢ to re,~oh'c the inlcrl~*ce shape in Ihe I:i~,. 12 problem hut was ;ldcqu;*le tier the I:il,, 11 ptohlcm. An

UflSOlvcd p~x)hlcm is m d¢l¢lnlil1¢,wilhotlt glid l'¢Iincnlc;11wh¢ll lhc interJit¢¢Sll;ip~i~;Jshccn f¢~olwd. The I¢~ull'+ presented above clearly nhow th;~l lhc nl¢lhod proposed ill Ihis p;Ip~:r ¢;tn lcsolvc Ct+lllpli¢illCt| jlll¢l|)lc¢,+ ill two splice dirtlCnsiorls. It is ~lraigllll',)rw;jld 10 cxlCnd llljs II]clh()d to

lhl¢c H'mc¢ di111¢llsi(Ins.Nolo th;]l I]1¢ ii11¢I|'ac¢bOUlld;il~ J~i confined wilhfll a cell. Two in(¢rl:Jce I~ouml;~ik-~, can ¢om¢ I£,g¢lhcr '~'ithin a +.'ell and cau,+¢ :~ pill'| 111IhC iill¢l!)g'¢ bollndaty hl ,£p;11';11¢ liom 11;¢ ICM, 'l'jfis is wllal h;q~pCncd in Fig,+, 1(i . n d 11, ']'h¢ level ,sol foimul:~li~l~

assulcs thal such ch;ulgcs ill1opolol,y ;~;¢h:~lldlcdwilhoul diflicully.I"b¢ inl¢ll)iccsCOllSid¢iT,'d in lhis ll;ipcrw¢i'¢ COlll;icldi~¢onlillgilics.|IIlh~ h=lui¢ we will ;d!~1111111ocxlcnd lhk mclh~d | 0 ~])¢;ll' I;ly¢lS,

AppcndlxA For 111c ,%;]k¢ of' o}rllpicl¢ll¢~,5 wc prc~¢nl ;1 ~,c(.'ond-oFdci i;xt~n~,JOll ~)l thc cla'~qc tirkt-old~.'J ~(.'h¢l'Kic o[ (~)ut';JM|. J~r~aC~()n ;~n {J }~¢C~ i T] ' Tiff,, is lhc m¢lhod ,,vhich i,,. tthcd |0 ~,01'¢C(11¢ Ic'vcI '~el equation i~k, -- I~

;~1

'~'~, = O.

(A,])

,lflfl

S,I". Dal i~ / 7~ ic L ',g,lor h~;~,'dmlic .~.¢~t='flll

As described Ilelow il is. a scalar, mln.o.n~sc~,ative ve [email protected] of Ihc conservative van l.eer-I Lmcock [33] MUSCL scheme. Consider a one space diincnsion version of equation ( A I ) , ilain¢ly fig fig -lq- + u . . ,'kv ..

O.

(A.2}

I.~ividc the real line into cells i.llld aplmlximale g by a linear flll!¢liOll till each cell, Th;It is G( x, t) = G( a,. t) + {x -.t, )S~ for ~ _ ,/. ~ . t < .~.'~, i/:,

(A.3)

Ilcrc ,'¢~ and G(x,, t) rcprcscllt aplmlxinmtions Io tile slope and average vah,e of g(x. t) ,.lrl cell x~, lcspecliv¢ly, G(x, l) is an ;ipproxinlalhln I0 g(.l', I), Ill Ih¢ folhv~'ing we Icl (IS delltqe G(. L, 1,,). T h e ~,ohlliOI1 lilt cell x ~11 1iln¢2 I,~, i :, is colnr)ulLd using ib¢ folllltll;i

c;;" ' / :

= ~;;' - Izy,"h,.S',.

(a.a)

I|¢re h I ~:.it' , I/~-'l'l I/2 ;lIid ilk):= ml/h F Ill 111111~ space dilncnsitlllS, eqll;Iti~lll (A,4) wttuld CI;lll;lill addilitmal ICI'IIY~ involvillg vt'h~cily ¢Olllponc'll~, lllld slopes ill (Hilt'l dircclions. ('ell inter fate v;llues arc computed using tile [OIIItU]:IC, • ~t, I//2

( , l ~, ~ ,, ~ (;;" ~ : "fl~l,I/~

4

(A.5)

',h,S~,

I;.I

I"

finally, lilt. ~,(lhlliOll ;11 ;ill1(,? I n , I iN k'¢llnr~tilcd tl:~illl~ Ih12 liffllltlla

~;;" I : (;;'.-- ~,U,"'

I/2

"*t ~ I/~, (~'n*#/1/211: 1(,,+,/: -

,

(A.6)

•~ ht3c (~..

[G/., till,

it' U > 0, ilU<0.

(A.7)

Apllrndlx II Ill thk appe.ndix w¢ dclJv,c the :q)p;oxl,;,4:]t¢ riCnl;llln r,ohHiofl Ihal v¢:~.~tl~,¢d to d¢lcl'lltltl¢ tile ~ll¢ll~th ~11 Ib¢ iIH¢II;ICe :~l)tl I0 t'~nlllule nUlllClical IILiXC',, lhtou~t!h the IloutIdaliC~, o~ cclF, cc~;~r;.ining inb:i|:]ces, 'l'hi,; appr,oxim=l¢ Riemann ~,~flulion i'~ trom a da',,s ~l appr~lximation~ 1hal wa,, lirq I~rcq'~t~,cd Ily I|arlcn, I.ax mid 'val~ l ¢¢r [|7]. Tller.¢ ;.;pl~tOximafion~ w.cl¢ u~cd lo dcliv¢ nunlclical tlu~ I,olmula¢ fill i?encl'al :,y',lcm', ol coll,ser,,,aliori laws by I);w!s [1/)], "HIe :~ppl~l~,,inlal¢ Ricmann s~luth)rR is dclivcd Ily a~,umin~.~ lhal an initial q,~, nlinuiiy in Ihe '.tlluthln bleaL~, up into k w;Jvc~,, ¢;~cll moving ;Jr a diflebcnl vch~cily aj. ,2 ..... k, Thi.~ ~llualilln i'~ iltuq!'aled in I:if., ILl, If we i~legt'al¢ Iht! t,on~¢lv;llion laws, equal,,~11 (3,1). over lh¢ Iccl;Jnkq~ ( - ~ti. i,/I) >: (I. ! 4 At). Ib¢ IO,0!| i~

IOg,)-/{u,

)-.,(l,~ -.,]

~ ~:{.~ - 1 4 ) ~. . . . ~ a d . , . . - , ~

,),

0~.~)

S.I:. D.ri~ t 7)ac!,ing/orIwperb#[i(" ,V,'sl'ms

467 u~ t

~l)_l

i

/ l~ I/l J

/ /:

~i,

u]I

Fi~!,I|.l. APpl(llilnal¢I(iclllalltlsolutionInod¢l.

/

/u;

ajat l:i~. B." Modclof;m isolated Iravcling tl~Jvc.

where

,; " (Ij, )~t -aj~! j,*~ar .(.~,1+aOdx.

(r).2)

Formula (ILl) defines a ilia" diJfcr¢lu'e .v)lilling (cf. [17]). If we itbtcgrat¢ (3.1) across an isolated wave Mth vci(:cily aj (sc¢ Fil,, B,2), the result is

wllcr¢

L' - f : , "",(": - ' ; ,),

0).3)

f" ~ ~ Jr

(n.4)

f(,(X, t)),.It.

EqlmliOI) (B.3) is ~) R:mkinc-llu~.oniol condJ|iol) which rcl;~Ics the jump in .~pace-avcra~cd v;))iahlc~ u 7 Io Ih¢ jump in fin)c,avcrai-tcd f[tlX¢~ f,'~ and Ihc wt)v¢ vclocily a r l| is (o b¢ ¢xl)cclcd I]):tt i)1 gcncl'al

fJ*~f('7)'

(B..~)

We nolc in pas%ing thai thi~ Iormulalion iJ)cludcs Ihc )u)pul~r r o c [27] mclhod, "I'<)scc Ihi~, wc tccall (hat Roe was ablc to conslrucl a mcan v:~luc inattix A(u=., uj~) !0r Ihc Eulcr cqua|ions whh ~n /deal gas cqualh!ll of ),lalc, Among olhcr Ihings, Ihis mallix sa(.i~ficd the. cottsclvalio)! c(mdifion

J(u,t) ' f O g . } - a ( u l , ~'r)(u),.- u)).

(IL(,)

S.F. Daris / 7h.'/;hlg for h~Tlerboh'c,i)'.lrcnl~

46,'~

'l'h¢ £)rmulati,m (B.I) billows ilnmedi:ltely from (11.61 with k equal to tile number or equ:llions in !!1e sys!,_'n~ provided ..... ~ (;t~ - i t ; _ I) :.'re lakcn ,o be the projeclions of ( % ~ - "L) on Ihe rigid eit'envectors of A(Ul., nil) and aj are lak'Jll Io b¢ the corresponding dgenvalucs. Ilarlcn, Lax and wm Leer [17] have shown that a firsl-i]rder melhod based ou Ihe approxiluat¢ ];',ien'lanrh s,.)ltltiol| (B.]) satisfies an cnlropy ¢olld[iion if "q is less than or equal to the minimuul wave speed ill the ex[icl Ricmanu solution and a.~ is [;re.'|tcr than or equal to lh¢ maximum wave :.,peed in lhc exact Iz.iem;uln soluliun. They a]s,a .,,llow lhat if only lwo waves are used in the appr¢~xinlation, lh¢ flu× diflcrerv..'c sp]illing takes lhc f,tlrm f ( U r ) - d ( u L ) =eL(u* uL) +ar(uj, ~ - u * ) . (I].71 Once lhe wave speeds a L and a~ arc specified, simple dosed form expressions can be obl~dned h)r u*, Delails are shown in Appendix C. This m,Jlh,c]d call Ix.: in|urpleled ;is a Godul)OV-lypc sch¢llle based or~ ~ll] approximate Ricm:mn slllulion wilh one inlermedi:~Je slate u *. The inler!r,e,!i~,_!e .,.!:~!¢ :!]:pr'.:::i:,-.a!es a space u'..,¢rag¢ ()f Ihc hllerlnt2dialu slIlles in the ¢xilcl I>,iemtlnrl solulk)n. The paper by I)avis [10] prupo~,cs a mmlber ~ff possible choices for lhc wave speeds el. and rill alld presenls llUnlCrieal r¢,~ults. "l'hcsc rcsulls show that Ihis sinlplc approximalion produces a¢¢uralu results in s¢cond-t~rder method,;. Therefore it is nul nccessaly Io use a more complicated illtttllt)d o]1 cell boundaries Ih,il do 111)1separate different malerials or trucked conliict dis¢ontinuilies. If the inlcrfae¢ being trucked lie.,i between lwo grid points, an averaged inlermedial¢ state would allow mixed malerial Io cr~ss the cell boundu~'. This would ¢reale a mixed malcrial ¢cll a i d deslroy lhe subcell resoluli(m sehenle. To avoid this we con~,lruct an approximate I;',ielllann solulloll Ihat contains ;,11 i,llerfac¢ and use it ollly ill inlerface cell boundaries. On such boundaries Ihe Ilux difference iplilling becomes

f(',3-/("D

= .,(i,i' - . L ) + V(,~

-'?)+",~('i~-uD,

lB.8)

llcr¢ V is tile speed of the malcrial interface and (u~ -u~') is llle jump ill lh¢ solulion across lh¢ malerial inlerfac¢. The qu:ullilie.~ are determined by Ihc meth(ld tlf chalacierislics (cf. [29]). The ~vsl of this :Jppcndix is u description of lids computation in (lie case where lh,." maleriuls satis[~, eilhcr file T~fit cquatioJ~ (3.8) or Ihe ideal gas equation ~f slal¢. The Tail etlu[llion of slille immcdi~,!dy implies Ihe relulion, (B.O)

( p + 1 ) / I , ~" = const.

This expres,.ion can be rearranged Io ,Jcfine Ihe COllM;lUl .~;~ I1y 2 I S~ ,~

In ¢'~ -- --

Int l) 4- 1 ),

(11.11))

where ,:',= l y ~ ( p + 1)//~.

lB.11)

For ihe ideal gas equalion of simile, sve assume Ihul Ihe pressure waves, Ihul is. th()se wilh speeds given by a L and a r are iscnlrqpi¢. This means thal we can replace Ihe c,crgy cqualion by lhe condili~m that lh¢ cnlrepy S~ c clined by 2 I

[~l~ Dat i~ / "l'l,wking fitr b.~l~clb~dic ~)~lcml

dfiq

is conManl, hi lh¢ above ~,~ is Ih¢ ralio of specific heats Ik~rall ideal gas aud c,~ = i'y~p/it is Ihe speed of ~;Ouild. The assunlplion thai S given Ivy either equatitm (B.10) or (I1.12) is constant tin C;lch side ol + Ihe inlerfi~ce implies Ihal the RiClllann invarianls, 2 R*= ~ t " +tt (B,13) and R=

2 "~._ ] - c - u .

(B.14I

are col~slant along d x / d l = u +c. m~d d x / d t = u - c respcclively, on cadt side of Ihe crfac_L w]t,2re c deJlote:; cilhcr ('~ or c~, its approprJale. Likewise, 'y dcllOlcs eitiler %, or Yt.a n d , is the velocity conlpOllCnt i1) the x-direction. Across the interface

v~ =

t,.,*

(i].ls)

u~ =tt~ = V

(11,1¢0

arid arc satisfied. These equations determine tile interlace veloeily 1/ and in|efface jump u~ -tt~. There arc a number of ways to solve these equali,tms. The ~Jpproach laken for this paper w;ts as follows: (I) Eliminme c~ between equations 113.13) and either (B,IO) or (11.12) as ~q~propriate. The result is an cqualian for p~ a', a funclion of itS. (2) l-lhnimJle c~ belween equalions (13.14) and either (B.IO)or (11.12);is appropriate 'HIe rcsull is an equation liar p~' as a funclion of u~. (3) Use equation (B.16) Io set V= u~ = u~. Suh'e cquali0n (B.15) in Ihe form

p~(I/) -p~(I/) = 0

0].17)

by Ncwlur~'s method for V. (4) Use equations (13.13) and (11.14) to delermin¢ c on each side of Ih¢: interl2~ce. ()blain P~ ~t'~ from eilher ¢qualilm (B.10) or (13.12). Then the ,,alucs of p on each side of Ihe intcthJc¢ can be delermined from Ihe definition tff c. In practice, the interfac,: julep is used to COIT¢Cl Ih¢ flux al the 1:.,3ur~dary between cells containing different malerials. Appendix C describes these flux correetioas, This same melhod of characteristics procedure is used to determine inlerfa¢¢ jumps which arc added 1o Ihe .'a~lution in a ecll when the level set fun¢lion changes sign.

Appendix C

7his appendix comains a description ,ff the numerical flux used fl~r Ihe calculalltms shown in (hi.'; paper, The basic scheme is des¢lihed in detail in [10]. 1tore we repeal Ihc delivalkm ~l' Ihe

|ltlmeri,:al flux that is used 11] IllosU p;:rls of the d(ml;lin Ihllt do 11['11contain inlcrfaccs. Then we derive corrections film are used at cell boundaries thlit scpar;ite different nuLterials. Away from any iutcrEtcc, the diffcrcucc in flux across ~1 cell bo.udary is split according to the formula (B.7). The eorrcsp{mding difference in suluticu) values is split according to

(.,(-.,.)

=

(.*

-.1)

~ (",t -"*)

(C.I)

These equations can be s(flvcd fl~r the waves (u* - u L) a,d (u u - . " ) .

The result is

("" - " l ) = [ / ( " p 0 - f ( " l ) - ' d u u - u L ) l / ( a l - a u ) . (u v. - ,1 * ) ~

[f(ul, ) -f(u,)

-

al.(ul¢ - " L)]/(rq(

- at. )'

(C.2) (C.3)

Using these results the uumeric~d flax I:(.l, ul~) is determined as follows: • If a~ >(L then

I"(.,. "u)=Y("I i. • else if at~ < 0, tllen

V(,t,...)-/(.,,).

• else P ( u l , . . ) - f ( u t ) + a d . " - .l. ) =f(um~ ) - au(Ut: - u*). If a cell bound;fly sup:mites rcgkms oonlai¢li;)g different materials, then Ih¢ flux difference is split :lccording to cquati,m (ITS). The lerm ( u [ - u ~ ) i s Ihe wave corresponding to the iutcrf;ic¢ and Iz is its velocity. Appendix B expl;~i.s how Ihey ;Ire computed. This interhlce is not permitted to chief a cell I,nundar~ since this wotJld erealc a mulli-material cell Thus Iwo i|umcrit'al fluxes I; I and Fit. ¢orrespondillg to condilions to Ibc left :rod righl of Ihe inlerface Icspeclively, are computed al such a cell bound;try The ,cw rlunlerieal fluxes can be eonlpuled as corr¢cljo,s Io Ihe non-interface fluxes. This is csp¢.:ially simple in Ihc two supersonic eases. Assume that aj. < V < alo Then: • If at.> 0,

I:l(",., "~) ~/"("J.. "~0, V:~(.,..j~) ~ I:(.,,

.,~)-- V(u~

-.~').

• If a ~ < (;,

V,(uL, .~) ~ I,'(.~.. "~0 - V(.~ .V,~(.l ,

.~ ),

..) "I:(uL, ..).

In the subsonic c;lse, we must account for the fact Ihal (u '~ - Ul.) ,/- (u~' - Ul.) and ( . . - u*) (u i~ - u~ ), To do this, we solve (ITS) along wilh the ~,plilting

(".

-", ) ~ (",~ - "~) + ("; - "~) + ("~ - "~)

alld use :qualions (C.2) arm (C,3) Io get

(";~-'q)~("'

. ~ + (V-a,t). •

.

(C.4)

S.F. Davi+ / TnrckingJor lop. Ihotic systcm~

471

and

(".-"h=(',,.-'*)

(c.6)

_+,.--~(u=. - (+,-(v-+~,.) *-ui%

T h e corresponding subsonic ilnmerica] flux formulae are: FI+(++I+ ' //I¢ ) = / ( U L ) + //l+(U? --IlL)

,:, . (~,~-v) =1 t.L, u . j - a ~ . ~ ( . ~

,

,

-., )

and Fu(Ut., lilt ) = f ( U r ) - att(u R -11~ ) =FOg.,

"") + ~"--,j,~--.~ (( "k '+- a-+, .() ~-t , ,4*

).

Ackno,Medgenmnl This work was supported by the NSWC Independent Research Fund. 1 would like to thank Dr. Stun U s h e r for suggesting that the level sel function be '~mooth and the referees for inspiring sonic additional thnughl on cite subcell resolution.

References II] V.A, Andr0nov, $.M. B;tkhrakh, E.E. Me~hkev. V.N. Mokhov, V.V. Nikiforov, A.V. Pevnilskii and A.I. Tol:dmtyakov, Tulbutunl mixing at Cent;tel sul[ace ~cceleratcd by shock waves, Savlel PpI)'~. JETP 44 (1976) 424-427. [2l W.|:, llallhaus and M, ]lull, Interaclinn Ilctwccn tile ocean suzl"accand underwater .,,phcrical blast waves, Phys. FI,idr 17 (1974) 1068. [3] P.K Burr :rod W,T. Ashurst, An inlclface schenle foi turbulent fl:~m¢ propagation+ ri~pl. SAND 82-8773 Sar~dia Nuti
472

S.F. Dells / "l"rack#lgf,~rh)7~erbolle.~y~tctns

[ 13] M,J. I:ritls |1lidJ.P. |]t,ris, The Laglarlgia. soltlliOllo| tr~lnsicnt problems in hydrodynamics using a [riangnhr rllc,,h, ]. Contpttt. Phy~'.31 (1979) 173-215. [14] J. Glirmn. E. [sa:lc:~on.D. Mard~csin and O./'4c||ly:m. Front Irackln,~ G)r hyperbolic s),ste~;, Adl. ht Al:pl. Magi. 2 (198I) 91-119, [15] A. Ilartcn. Iligh rcsdntion .',chctncs fnr hypcrll()liL:con,,err;ilion laws, J. Compul. Phks. 49 (1983) 357-393. [16] A. |larlen. ENO schclncs ,.vithsubcdl resnlutk)n, Z G:mlmt. Phyr. 83 (1989) 148-184. [17] A. |[atlcn, P.D. Lax and B. van Leer, Oil IJpMrC;llndifflzrcncing and Godun(w-lyf~¢ sch~nleg fill" hyperhoIic conser','alh~n laws, ShIM Rrr. 25 (1983) 35-61. [18] C.W. Ilirl and B,D. Niehnls, Vohmle iff fluid (VOF) i~lclla)d h:r the (lyt1;lrllicso:' free bt~undarics. ,L Comp, t. Ph)~. 39 (1981) 2111-225. [I9] J.M. llym:m. Nunleficnl methods filr Iracking interface,;, in: A.I,LBishop, L,J. C:,,!!pbdl and P,J, Channcll, ads., Fronls, htterfucc.~ and Path't~ls (Elsevier. New York. 1984). (2(I] K,M. Li and M. 11o11,Nurneri,:al .~ohlliollSI(I V,ralcr wnves gcll¢lalqd by Sllallt}wunderwater explosi(irls, Ph}w, Fluidt 24 (1981) ~116. [21] W. Muldcr, S, Oshcr and LA. Selhi:m. Computing inlcrfaec motion in compressible gas dynamics, CAM Rt:pt. 9(1-03, Uni'.'cr:,hyo| California, (a.. An~zclcs,CA (19911). [22] ll.D. Nit-ht~kand C.W, It~rl, Method~.filr calculating mulli-dimerlsional, tran,,ient frcc surface flllwsp',l~,lIlodics, Rcpl. LA-UR-75-1932 Los Alanlos Scicnlific [.:d~nr;ikllT.ta)s AhmloS.NM (1975). [23] W.F, Nob and P. Wot~d¢,ard. SI.IC (simple line inlerfa¢.: ¢~lculali(m),hi: P.J, Zandbcrgcn~ cd,. I'rom,dings of l],. Fif,'~,L;;cn;.thJttal f_)mfen.m'e on Numcfil'a! ~h'thod.~ hi Fluid D_vnamks ~Springer, Bed in, 1976). [24[ S. Oshcr i~nd J.A. St2lhi:ln, Frtlnls pl.p:lgaling '.'.ilh curvatlJrc-dt~pcfldenl speed: aJgorilhnls based nn ]tamillon-Jacobi filrmukttions, .L Cm~qnlt. I'llys. 79 (19~8) 12-49. 125] E,A. Overman :rod N.J. Zabusky. I-volulion ;rod mcrgt:r of isolated '.(:flex struclures. Repl. ICMA-81-3I ]rl~,lilul~.' J'OrC()mpulaliollal niHh¢lll;llics and Apl)Iicatit~ (1981). [20] |~.,D. RJfhJlJl]ft2[, Taylor in~.labiiJlyill sl)oek acccletati,.m i)f con,ptcssible Jluids. Comnt. Pllre AppL Math. 13 (19(,~)) 297-319, [271 I'.1.. Rot:. Apprt*ximalc Riem~nn ~k'cr... p~ranlclcr vectors,and tliffcrcac,.:s,:hcmcs. 3. Compt:t. Phys. 43 (1981) 357-372. [2~] P.L. Roe. Discontinuous solutinn,, to, hyperbolic ~:yslem'.under opt:ralor :,plilling. RcpI, No, 87-64 If'AS|', |tampl()n, VA (1987). [.'.9] G. Rudinger. II~tv Diagmmx for Non,toady Fhm' ht Dn~'ts(Dover, New York, 1969). [311[ C.-W. Shu and S. ()xl~er, Elficien[ implemtrntalio!l of,,:ssentially nan-tlscillalnq slmck ¢apluring schemes IF, J. Ctmlput. Ph)w, g3 (1989) 32-78, 131] G.A. Stsd, Nunwrica! 3fi'tlu;dv in Fhdd Ihnamirs (Cimlbridge Univcrsily ]kcss, Cambridge, England, 19i~5). [32[ G, Sir;rag. On Ih¢ c,.lnslruclion and comparison of difference schemes, SIAM J. Numer. Anal. 5 (1968) 5[)6~517. J33] G.D, :an Albada, B.','an Leer and W.W. Roberts A c~mpara t',: s dy fc n pu eli naI nlcthods in el ~mi,:ra~ dynamic:,, A~tr~mom. Astrophys. l(l~ (1932) "I6-84. 134] NN. YancnLo, 771eMctl~M ofFracthmal Stcpl (Springer, Berlin, 1971). [35] IF. Y;~nt~, An ailifit:iM¢~mprcssion method for ENO :,chcmlts:The slope modilic:~linnmclhtld, J. Comput. l'hy.~. 89 119911)125-lggL [36] I).L. Youn,~,,, Numclic;d simul:ltion (ff tnrbuicnl mixing I,y I{aylcigh-Taylor inslabilily, Ph).~. 121)(1984) 32-44.