## A high-resolution mapped grid algorithm for compressible multiphase ﬂow problems

### K.-M. Shyue

^{*}

Department of Mathematics, National Taiwan University, Taipei 106, Taiwan

a r t i c l e i n f o

Article history:

Received 3 May 2010

Received in revised form 16 July 2010 Accepted 8 August 2010

Available online 18 August 2010

Keywords:

Compressible multiphase ﬂow Fluid-mixture model Mapped grids

Wave-propagation method Stiffened gas equation of state

a b s t r a c t

We describe a simple mapped-grid approach for the efﬁcient numerical simulation of com- pressible multiphase ﬂow in general multi-dimensional geometries. The algorithm uses a curvilinear coordinate formulation of the equations that is derived for the Euler equations with the stiffened gas equation of state to ensure the correct ﬂuid mixing when approxi- mating the equations numerically with material interfaces. Ac-based and aa-based model have been described that is an easy extension of the Cartesian coordinates counterpart devised previously by the author[30]. A standard high-resolution mapped grid method in wave-propagation form is employed to solve the proposed multiphase models, giving the natural generalization of the previous one from single-phase to multiphase ﬂow prob- lems. We validate our algorithm by performing numerical tests in two and three dimen- sions that show second order accurate results for smooth ﬂow problems and also free of spurious oscillations in the pressure for problems with interfaces. This includes also some tests where our quadrilateral-grid results in two dimensions are in direct comparisons with those obtained using a wave-propagation based Cartesian grid embedded boundary method.

Ó 2010 Elsevier Inc. All rights reserved.

1. Introduction

Our goal is to describe a simple mapped-grid approach for efﬁcient numerical resolution of compressible multiphase ﬂow in general multi-dimensional geometries. As a ﬁrst endeavor towards the method development, we are concerned with a simpliﬁed model problem, where the ﬂow regime of interest is assumed to be homogeneous with no jumps in the pressure and velocity (the normal component of it) across the material interface separating two regions of different ﬂuid components within a spatial domain. In this problem, the physical effects such as the viscosity, surface tension, and heat conduction are assumed to be small, and hence can be ignored. With that, we use an Eulerian viewpoint of the governing equations that the principal motion of each ﬂuid component in a Cartesian coordinates can be written as

@

@t

### q q

uiE 0 B@

1
CA þX^{N}^{d}

j¼1

@

@xj

### q

uj### q

uiujþ pdijEujþ puj

0 B@

1

CA ¼ 0 ð1Þ

for i ¼ 1; 2; . . . ; Nd. Here Nddenotes the number of spatial dimensions. The quantities

### q

, uj, p, E, and dijare the density, particle velocity in the xj-direction, pressure, total energy, and the Kronecker delta, respectively.0021-9991/$ - see front matter Ó 2010 Elsevier Inc. All rights reserved.

doi:10.1016/j.jcp.2010.08.010

*Tel.: +886 2 3366 2866; fax: +886 2 2391 4439.

E-mail address:shyue@math.ntu.edu.tw

Contents lists available atScienceDirect

## Journal of Computational Physics

j o u r n a l h o m e p a g e : w w w . e l s e v i e r . c o m / l o c a t e / j c p

To close the model, for simplicity, the constitutive law for each ﬂuid phase is assumed to satisfy a linearized Mie-Grün- eisen (i.e., the linearly density-dependent stiffened gas) equation of state of the form

pð

### q

;eÞ ¼ ð### c

1Þ### q

e þ ð### q

### q

_{0}ÞB ð2Þ

for approximating materials including compressible liquids and solids (cf.[14,25]). Here e represents the speciﬁc internal energy. The quantities

### c

,### q

0, and B are the ratio of speciﬁc heats (### c

> 1), the reference values of density, and speed of sound squared, respectively. We have E ¼### q

e þPNdj¼1

### q

u^{2}

_{j}=2 as usual.

In this work, we want to generalize a state-of-the-art shock-capturing method that was devised originally for single-phase ﬂows on mapped grids to the case of a multiphase ﬂow. It is well known that the principal problem in the usual extension is the occurrence of spurious pressure oscillations when two or more ﬂuid components are present in a grid cell (cf.[5]and references therein). Here the algorithm uses a curvilinear coordinate formulation of a ﬂuid-mixture model that is composed of the Euler equations of gas dynamics for the basic conserved variables and an additional set of effective equations for the problem-depen- dent material quantities. In this approach, as in its Cartesian coordinates counterpart (cf.[30,31]), the latter equations are de- rived to ensure the correct ﬂuid mixing when approximating the equations numerically with material interfaces, see Section2.

With the proposed model equations, accurate results can be obtained on a mapped grid using a standard method, such as the high-resolution wave- propagation algorithm for a single-phase ﬂow (cf.[8,9,21]), see Section4for numerical examples.

There are quite a few other numerical approaches available in the literature for approximating compressible multiphase problems over a multi-dimensional domain with complex geometries. Some representative ones are the overlapping grid method[5], the unstructured grid methods[2,11,43], and the Cartesian grid embedded boundary method[35].

An advantage of the mapped-grid approach described here is that extension of the method from two to three dimensions can be done in a straightforward manner for simple geometries such as cylinders, spheres, and their variants[9]. This is in contrast with the extension of an unstructured or a Cartesian grid method, where it requires a signiﬁcant algorithmic and programming effort to realize each of the methods that are designed of general purpose. In addition to that, if we want to use a front tracking method to improve numerical resolution of shock waves and interfaces, with complex geometries involved, it would be relatively easier to apply the method on a body-ﬁtted mapped grid than on a ﬁxed Cartesian grid, see[17,37]for an example. Furthermore, it is also easy to combine the method with a class of moving mesh techniques (cf.[38,40]) for efﬁcient solution adaptation.

It should be mentioned that the methodology we have given here is by no means limited to the current case with the stiffened gas equation of state. Extension of the method to problems involving more complicated equations of state can be made by considering, for instance, either a

### c

-based model of the author[32]or a### a

-based model of Allaire et al.[3](see Section2.2), and proceeding with the idea described in this paper. Without going into the details for that, our goal is to establish the basic solution strategy and validate its use via some sample numerical experimentations; this is a necessary step for our further development of the method towards more complicated problems of fundamental importance (cf.[15,18,27,28]and references therein).

The format of this paper is as follows: in Section2, we describe our mathematical model for a simpliﬁed homogeneous multiphase ﬂow in curvilinear coordinates. In Section3, we review brieﬂy the wave-propagation method on mapped grids.

Numerical results of some sample test problems in two and three dimensions are presented in Section4.

2. Mathematical models in curvilinear coordinates

The basic governing equations in our mapped grid algorithm consist of two parts. We use the Euler equations in a curvilinear coordinate as a model system for the motion of the ﬂuid mixtures of the conserved variables in a multiphase grid cell. With that, from the mass and energy conservations, we derive a set of effective equations for the problem-dependent material quantities in those cells, see below, that can be used directly to the determination of the pressure from the equation of state. Combining this two set of the equations together with the equation of state constitutes a complete mathematical model that is fundamen- tal in our mapped grid algorithm for numerical approximation of multiphase ﬂow problems with complex geometries.

To ﬁnd out the aforementioned equations in a three-dimensional Nd= 3 curvilinear coordinate system, for example, we introduce a coordinate mapping from the physical domain (x1, x2, x3) to the computational domain (n1, n2, n3) via the relations

dx1¼ a1dn1þ a2dn2þ a3dn3; dx2¼ b1dn1þ b2dn2þ b3dn3; dx3¼ c1dn1þ c2dn2þ c3dn3;

ð3Þ

where ai, bi, cifor i ¼ 1; 2; 3 are the metric terms of the mapping. Then under this mapping, the Euler equation(1)can be transformed into the new coordinate system as

@

### q

@t þ1 J

X^{N}^{d}

j¼1

@

@nj

ð

### q

UjÞ ¼ 0;@

@tð

### q

uiÞ þ1 JX^{N}^{d}

j¼1

@

@nj

ð

### q

uiUjþ pJjiÞ ¼ 0; i ¼ 1; 2; . . . ; Nd;@E

@tþ1 J

X^{N}^{d}

j¼1

@

@nj

ðEUjþ pUjÞ ¼ 0;

ð4Þ

where Uj¼PN_{d}

i¼1uiJ_{ji}is the contravariant velocity in the nj-direction for j = 1, 2, . . . , Nd. Here the quantities Jijfor i, j = 1, 2, 3 are
as a consequence of the coordinate transformation that satisﬁes the following expressions:

J_{11} J_{12} J_{13}
J_{21} J_{22} J_{23}
J_{31} J_{32} J_{33}
0

B@

1 CA ¼

b2c3 b3c2 a3c2 a2c3 a2b3 a3b2

b3c1 b1c3 a1c3 a3c1 a3b1 a1b3

b1c2 b2c1 a2c1 a1c2 a1b2 a2b1

0 B@

1

CA ð5Þ

and the quantity J = detj@(x_{1}, x_{2}, x_{3})/@(n_{1}, n_{2}, n_{3})j is the Jacobian of the mapping which can be computed by

J ¼X^{3}

i¼1

aiJ_{1i}¼X^{3}

i¼1

biJ_{2i}¼X^{3}

i¼1

ciJ_{3i}: ð6Þ

Note that during the initialization step, all the coordinate transformation variables such as ai, bi, ci, J1i, J2i, J3ifor i = 1, 2, 3, and J would be determined and remained ﬁxed at all time when a mapped grid is constructed by a chosen numerical grid gener- ator (cf.[9,39]).

It is easy to see that(3)would be a two-dimensional coordinate mapping from (x1, x2) to (n1, n2) for any spatial location x3

in the physical domain, if we have a simpliﬁed data set where the quantities a3, b3, c1, and c2are all zero, and c3is equal to one. In this instance, if we set Nd= 2 in(4)with the coordinate transformation variables deﬁned as in(5) and (6), we would have the same Euler equations in a two-dimensional curvilinear coordinate when a mapping of the form

dx1¼ a1dn1þ a2dn2;

dx2¼ b1dn1þ b2dn2; ð7Þ

is used in the derivation (cf.[4,7,16,42]). Thus, without causing any confusion, we may simply use the symbol Ndas in the Cartesian case, see(1), to represent the number of spatial dimension in the curvilinear coordinate formulation of equations.

2.1.

### c

-based model equationsTo derive the effective equations for the mixture of material quantities in curvilinear coordinates, one approach is to start with an interface-only problem (cf.[30,32–34]) where both the pressure and each phase of the particle velocities are con- stant in the domain, while the other variables such as the density and the material quantities are having jumps across some interfaces. Then, from(4), it is easy to obtain an equation for the time-dependent behavior of the total internal energy as

@

@tð

### q

eÞ þ1 JX^{N}^{d}

j¼1

Uj

@

@nj

ð

### q

eÞ ¼ 0:Now, by inserting(2)into the above equation, we ﬁnd an alternative form:

@

@t p

### c

1### q

### q

_{0}

### c

1 B

þ1 J

X^{N}^{d}

j¼1

Uj

@

@nj

p

### c

1### q

### q

_{0}

### c

1 B

¼ 0; ð8Þ

that is in relation to not only the pressure, but also the density and the material quantities

### c

,### q

0, and B.In our algorithm, to maintain the pressure in equilibrium as it should be for this interface-only problem and also to deter- mine all the three material quantities in (2), we split (8) into the following three equations for the ﬂuid mixture of 1=ð

### c

1Þ;### q

B=ð### c

1Þ, and### q

_{0}B=ð

### c

1Þ as@

@t 1

### c

1

þ1 J

X^{N}^{d}

j¼1

Uj

@

@nj

1

### c

1

¼ 0; ð9Þ

@

@t

### q

B### c

1

þ1 J

X^{N}^{d}

j¼1

Uj

@

@nj

### q

B### c

1

¼ 0; ð10Þ

@

@t

### q

_{0}B

### c

1

þ1 J

X^{N}^{d}

j¼1

Uj

@

@nj

### q

_{0}B

### c

1

¼ 0; ð11Þ

respectively. As before (cf.[31–34]), since in practice we are interested in shock wave problems as well, we should take the equations, i.e.,(9)–(11), in a form so that all the three material quantities remain unchanged across both shocks and rare- faction waves. In this regard, it is easy to see that with 1/(

### c

1) and### q

_{0}B=ð

### c

1Þ governed in turn by(9) and (11), there is no problem to do so (cf.[1,30]). For### q

B=ð### c

1Þ, however, due to the dependence of the density term, it turns out that, in a time when such a situation occurs, for consistent with the mass conservation law of the ﬂuid mixture in(4), the prim- itive form of(10)should be modiﬁed by@

@t

### q

B### c

1

þ1 J

X^{N}^{d}

j¼1

@

@nj

### q

B### c

1Uj

¼ 0; ð12Þ

so that the mass-conserving property of the solution in the single-phase region can be acquired also.

In summary, combining the Euler equation(4)and the set of effective equations:(9), (11), and (12), yields a so-called

### c

- based model system as@

### q

@t þ1 J

X^{N}^{d}

j¼1

@

@nj

ð

### q

UjÞ ¼ 0;@

@tð

### q

uiÞ þ1 JX^{N}^{d}

j¼1

@

@nj

ð

### q

uiUjþ pJjiÞ ¼ 0; i ¼ 1; 2; . . . ; Nd;@E

@tþ1 J

X^{N}^{d}

j¼1

@

@nj

ðEUjþ pUjÞ ¼ 0;

@

@t 1

### c

1

þ1 J

X^{N}^{d}

j¼1

Uj

@

@nj

1

### c

1

¼ 0;

@

@t

### q

0B### c

1

þ1 J

X^{N}^{d}

j¼1

Uj

@

@nj

### q

0B### c

1

¼ 0;

@

@t

### q

B### c

1

þ1 J

X^{N}^{d}

j¼1

@

@nj

### q

B### c

1Uj

¼ 0;

ð13Þ

this gives us Nd+ 5 equations to be solved in total that is nicely independent of the number of ﬂuid phases involved in the problem. With a system expressed in this way, there is no problem to compute all the state variables of interests, including the pressure from the equation of state

p ¼ E
PN_{d}

i¼1ð

### q

uiÞ^{2}

2

### q

^{þ}

### q

B### c

1

### q

0B### c

1

" #,

1

### c

1

: ð14Þ

For the ease of the latter discussion, it is useful to write(13)into a more compact expression by

@q

@tþ1 J

X^{N}^{d}

j¼1

@

@nj

fjðqÞ þ BjðqÞ@q

@nj

¼ 0; ð15Þ

with

q ¼

### q

;### q

u1; . . . ;### q

uN_{d};E; 1

### c

1;### q

0B### c

1;### q

B### c

1T

;

fj¼

### q

Uj;### q

u1Ujþ pJj1; . . . ;### q

uN_{d}Ujþ pJj;Nd;EUjþ pUj;0; 0;

### q

B### c

1UjT

; Bj¼ diagð0; 0; . . . ; 0; 0; Uj;Uj;0Þ;

ð16Þ

for j ¼ 1; 2; . . . ; Nd. Note that in the Cartesian coordinates case where the coordinate mapping quantities a1, b2, c3are all equal to one, while the remaining ones are all zeros, Eq.(15)reduces to

@q

@tþX^{N}^{d}

j¼1

@

@xj

fjðqÞ þ BjðqÞ@q

@xj

¼ 0; ð17Þ

with fjand Bjdeﬁned in turn by

fj¼

### q

uj;### q

u1ujþ pd1j; . . . ;### q

uN_{d}ujþ pd3j;Eujþ puj;0; 0;

### q

B### c

1ujT

; Bj¼ diagð0; 0; . . . ; 0; 0; uj;uj;0Þ:

ð18Þ

Then it is easy to check that fjand Bjare related to fjand Bjvia

fj¼X^{N}^{d}

i¼1

fiJ_{ji} and Bj¼X^{N}^{d}

i¼1

BiJ_{ji};

respectively.

With these notations, by assuming the proper smoothness of the solutions, the quasi-linear form of our model(15)can be written as

@q

@tþ1 J

X^{N}^{d}

j¼1

ðAjðqÞ þ BjðqÞÞ@q

@nj

¼ 0; ð19Þ

where Aj¼ @fj=@q ¼PN_{d}

i¼1AiJ_{ji}is the Jacobian matrix of fjwith Ai¼ @fi=@q for i ¼ 1; 2; . . . ; Nd. If we assume further that the
thermodynamic description of the materials of interest is limited by the stability requirement, it is a straightforward matter
to show that any linear combination of the matrices Aiþ Bifor i ¼ 1; 2; . . . ; Ndis diagonalizable with real eigenvalues and a
complete set of linearly independent right eigenvectors (cf.[33]). Hence, we may conclude that our multiphase model is
hyperbolic. Regarding discontinuous solutions of the system, such as shock waves or contact discontinuities, we ﬁnd the
usual form of the Rankine–Hugoniot jump conditions across the waves (cf.[13]).

2.2.

### a

-based model equationsBefore proceeding further, it should be mentioned that to deﬁne the initial ﬂuid mixtures 1=ð

### c

1Þ;### q

B=ð### c

1Þ, and### q

_{0}B=ð

### c

1Þ in a grid cell that contains MfP1 different ﬂuid phases where each of them occupies a distinct region with a volume-fraction function### a

i2 [0, 1] in relation to it for i ¼ 1; 2; . . . ; Mf;PMfi¼1

### a

i¼ 1, we use the equation of state(2)of the form### q

e ¼X^{M}

^{f}

i¼1

### a

i### q

iei¼X^{M}

^{f}

i¼1

### a

ipi

### c

_{i}1

### q

_{i}

### q

_{0;i}

### c

_{i}1 Bi

¼ p

### c

1### q

### q

_{0}

### c

1 B:Here the subscript ‘‘i” denotes the state variable of ﬂuid phase i. By taking a similar approach as employed in Section2.1for the derivation of the

### c

-based effective equations it comes out readily a splitting of the above expression into the relations:1

### c

1¼X^{M}

^{f}

i¼1

### a

i### c

_{i}1;

### q

B### c

1¼X^{M}

^{f}

i¼1

### a

i### q

_{i}Bi

### c

_{i}1;

### q

_{0}B

### c

1¼X^{M}

^{f}

i¼1

### a

i### q

0;iBi### c

_{i}1; ð20Þ

where in the process of splitting the terms we have imposed the condition

p

### c

1¼X^{M}

^{f}

i¼1

### a

ipi### c

_{i}1: ð21Þ

Clearly when each of the partial pressures is in an equilibrium state within a grid cell, in conjunction with the ﬁrst part of (20), the pressure p obtained from(21)would remain in the same equilibrium as well, i.e., p = pifor i ¼ 1; 2; . . . ; Mf.

Now if the above volume-fraction notion of the states 1=ð

### c

1Þ;### q

B=ð### c

1Þ, and### q

_{0}B=ð

### c

1Þ are being employed in the### c

- based effective equations together with the usual deﬁnition of the mixture density### q

¼P^{M}f

i¼1

### a

i### q

_{i}, we are able to rewrite them straightforwardly into a componentwise form as

@

@t

### a

i### c

_{i}1

þ1 J

X^{N}^{d}

j¼1

Uj

@

@nj

### a

i### c

_{i}1

¼ 0; ð22Þ

@

@t

### a

i### q

iBi### c

_{i}1

þ1 J

X^{N}^{d}

j¼1

@

@nj

### a

i### q

iBi### c

_{i}1Uj

¼ 0; ð23Þ

@

@t

### a

i### q

_{0;i}Bi

### c

i 1

þ1 J

X^{N}^{d}

j¼1

Uj

@

@nj

### a

i### q

_{0;i}Bi

### c

i 1

¼ 0; ð24Þ

i ¼ 1; 2; . . . ; Mf. Then based on the fact that all the material quantities

### c

_{i};Bi, and

### q

0,iwill be kept as a constant in each phase of the domain at all time, from(22) or (24), it is easy to ﬁnd the transport equation for the volume-fraction### a

ias@

### a

i@t þ1 J

X^{N}^{d}

j¼1

Uj

@

### a

i@nj

¼ 0; ð25Þ

whereas, from(23), we ﬁnd the conservation law for the phasic density

### a

i### q

ias@

@tð

### a

i### q

iÞ þ1 JX^{N}^{d}

j¼1

@

@nj

ð

### a

i### q

iUjÞ ¼ 0: ð26ÞIt is apparent that, if the solutions of

### a

iand### a

i### q

iare known from the equations for i ¼ 1; 2; . . . ; Mf, we may, therefore, com- pute 1/(### c

1),### q

B=ð### c

1Þ, and### q

_{0}B=ð

### c

1Þ directly according to(20). Thus, instead of using the### c

-based effective equations, it is a viable alternate to use the### a

-based equations:(25) and (26), for the motion of the mixture of the material quantities of the problem.To sum up, combining(25) and (26)with the momentum and energy equations in(4)yields a

### a

-based (or called volume- fraction) model system that can be written as@

@tð

### a

i### q

_{i}Þ þ1 J

X^{N}^{d}

j¼1

@

@nj

ð

### a

i### q

_{i}UjÞ ¼ 0; i ¼ 1; 2; . . . ; Mf;

@

@tð

### q

uiÞ þ1 JX^{N}^{d}

j¼1

@

@nj

### q

uiUjþ pJji

¼ 0; i ¼ 1; 2; . . . ; Nd;

@E

@tþ1 J

X^{N}^{d}

j¼1

@

@nj

ðEUjþ pUjÞ ¼ 0;

@

### a

i@t þ1 J

X^{N}^{d}

j¼1

Uj

@

### a

i@nj

¼ 0; i ¼ 1; 2; . . . ; Mf 1;

ð27Þ

this gives us totally 2Mf+ Ndequations to be solved. Here analogously to(14)the pressure can be determined from the equa- tion of state

p ¼ E
PN_{d}

i¼1ð

### q

uiÞ^{2}

2

### q

^{þ}

X^{M}^{f}

i¼1

### a

i### q

iBi### c

_{i}1X

^{M}

^{f}

i¼1

### a

i### q

_{0;i}Bi

### c

_{i}1

" #,

X^{M}^{f}

i¼1

### a

i### c

_{i}1; where we have assumed

### a

Mf ¼ 1 PMf1i¼1

### a

i.It is clear that(27)can be written of the form(15)in which we have q, fj, and Bjdeﬁned by

q ¼

### a

1### q

_{1}; . . . ;

### a

Mf### q

_{M}

f;

### q

u1; . . . ;### q

uNd;E;### a

1; . . . ;### a

Mf1T

;

fj¼

### a

1### q

1Uj; . . . ;### a

Mf### q

M_{f}Uj;

### q

u1Ujþ pJj1; . . . ;### q

uNdUjþ pJj;N_{d};EUjþ pUj;0; . . . ; 0T

; Bj¼ diagð0; . . . ; 0; . . . ; 0; Uj; . . . ;UjÞ:

In a similar manner, in the Cartesian coordinates case where the system is of the form(17), we have fjand Bjdeﬁned by

fj¼

### a

1### q

_{1}uj; . . . ;

### a

M_{f}

### q

_{M}

_{f}uj;

### q

u1ujþ pdj1; . . . ;### q

uN_{d}ujþ pdj;N

_{d};Eujþ puj;0; . . . ; 0T

; Bj¼ diagð0; . . . ; 0; . . . ; 0; uj; . . . ;ujÞ:

We note that since the derivation of the

### a

-based model follows closely to the### c

-based model, it can be shown that this model is hyperbolic also (cf.[3]for the Cartesian coordinates case) and is as effective as the### c

-based model for multiphase ﬂow problems with the stiffened gas equation of state. But for problems with MfP3, the### c

-based model is a prefer one to use, because the basic equations for the model stay as Nd+ 5, see(13), irrespective of the number of ﬂuid phases involved in the problem.2.3. Include source terms

To end this section, we comment that if x_{1}is the axisymmetric direction, an axisymmetric version of our multiphase mod-
els in two dimensions can be written as

@q

@tþ1 J

X^{2}

j¼1

@

@nj

fjðqÞ þ BjðqÞ@q

@nj

¼ wðqÞ; ð28Þ

wherewis the source term derived directly from the geometric simpliﬁcation. That is, we ﬁnd

w¼ 1 x1

### q

u1;### q

u^{2}

_{1};

### q

u1u2;Eu1þ pu1;0; 0;### q

B### c

1u1T

; ð29Þ

when the

### c

-based model is considered, and havew¼ 1 x1

### a

1### q

_{1}u1; . . . ;

### a

M_{f}

### q

_{M}

_{f}uM

_{f};

### q

u^{2}

_{1};

### q

u1u2;Eu1þ pu1;0; . . . ; 0T

;

when the

### a

-based model is considered. In addition to that, if gravity is the only body force in the problem formulation, in the### c

-based model, for example, we may add in the following source term:w¼ ð0; 0;

### q

g;### q

gu2;0; 0; 0Þ^{T}

as well. Here g denotes the gravitational constant. As to the other source terms such as the one arise from the surface tension force at the interface, we may use a continuum surface force model of Brackbill et al.[6]for that, see the work done by Peri-

gaud and Saurel[26]and the references therein for more details. Since it is beyond the scope of this paper, we will not dis- cuss this further.

3. Numerical approximation on mapped grids

We use a state-of-the-art ﬁnite volume method in wave-propagation form (cf.[9,21]) for the numerical discretization of our multiphase ﬂow models (without the source terms) on mapped grids. The method is based on solving one-dimensional Riemann problems at each cell edge, and the waves (i.e., discontinuities moving at constant speeds) arising from the Rie- mann problem are employed to update the cell averages in the cells neighboring each edge.

To review the basic idea of the method, we consider the two-dimensional quadrilateral-grid case as illustrated inFig. 1, for example. In a ﬁnite volume method, the approximate value of the cell average of the solution q over the (i, j) th grid cell at a time tncan be written as

Q^{n}_{ij} 1
MðCijÞ

Z

Cij

qðx1;x2;tnÞdx1dx2¼ 1

### j

ijDn1Dn2Z bCij

qðn1;n2;tnÞdn1dn2;

where Cijand bCijdenote the regions occupied by the grid cell (i, j) in the physical and computational domains, respectively, and MðCijÞ ¼

### j

ijDn_{1}Dn

_{2}is the measure (area) of Cij. Here

### j

ij¼ JðCijÞ is the Jacobian of the mapping of this cell, andDnkis the mesh size in the nk-direction of the computational domain for k = 1, 2. The time step from the current time tnto the next tn+1is denoted byDt.

In this setup, a fully discrete version of the wave-propagation method for the Eq.(15)is a Godunov-type scheme on a quadrilateral grid that can be written as

Q^{nþ1}_{ij} ¼ Q^{n}ij 1

### j

ijDt Dn1

A^{þ}1DQ_{i1=2;j}þ A^{}1DQ_{iþ1=2;j}

1

### j

ijDt Dn2

A^{þ}2DQ_{i;j1=2}þ A^{}2DQ_{i;jþ1=2}

: ð30Þ

Here A^{þ}_{1}DQi1=2;j;A^{}_{1}DQiþ1=2;j;A^{þ}_{2}DQi;j1=2, and A^{}_{2}DQi;jþ1=2 are the right-, left-, up-, and down-moving ﬂuctuations, respec-
tively, that are entering into the grid cell. To determine these ﬂuctuations, we need to solve the one-dimensional Riemann
problems normal to the cell edges.

3.1. Computing ﬂuctuations

Considering the ﬂuctuations A^{}_{1}DQi1=2;jarising from the edge (i 1/2, j) between cells (i 1, j) and (i, j), for example. This
amounts to solve a Cauchy problem in the n1-direction that consists of

@q

@tþ1 J

@f1ðqÞ

@n1

þ B1ðqÞ1 J

@q

@n1

¼ 0;

as for the equations and the piecewise constant data

qðn1;n2;tnÞ ¼ Q^{n}_{i1;j}; if n1<ðn1Þ_{i1=2};
Q^{n}_{ij}; if n1>ðn1Þ_{i1=2};
(

as for the initial condition at a time tn. As an example, we describe next the case for the

### c

-based model(13)in more details, and that the case for the### a

-based model(27)can be constructed in an analogous manner.Let ~ni1=2;j¼ ð^b2;^a2Þ_{i1=2;j}and~ti1=2;j¼ ð^a2; ^b2Þ_{i1=2;j}be the unit normal and tangential vectors to the cell edge (i 1/2, j) in
the physical grid, where ^ai¼ ai=Si and ^bi¼ bi=Si are the scaled version of the metric elements ai and bi in (7) with
Si¼

ﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
a^{2}_{i} þ b^{2}_{i}
q

for i = 1, 2. Then to compute A^{}_{1}DQ_{i1=2;j}, as in[9,21], we perform the following steps:

(1) Transform the data Q^{n}_{i1;j}and Q^{n}_{i;j}into the new data QLand QRvia
QL¼ Ri1=2;jQ^{n}_{i1;j}; QR¼ Ri1=2;jQ^{n}_{i;j}:

Fig. 1. A sample grid system in our two-dimensional numerical method on a quadrilateral grid. The numerical solution on the rectangular grid cell bCijin the computational domain gives distinctively the result on the mapped quadrilateral grid cell Cijin the physical domain for all the grid cell (i, j).

Here Ri1=2;jis a rotation matrix deﬁned by

Ri1=2;j¼

1 0 0 0

0 ð^b2Þ_{i1=2;j} ð^a2Þ_{i1=2;j} 0
0 ð^a2Þ_{i1=2;j} ð^b2Þ_{i1=2;j} 0

0 0 0 I

0 BB B@

1 CC CA

with I in it as being a 4 4 identity matrix. Clearly this rotation matrix rotates the velocity components of Q into components normal and tangential to the cell edge, and leaves the remaining components unchanged.

(2) Solve Riemann problem in the ‘‘x1” direction for

@q

@tþ @

@x1

f1ðqÞ þ B1ðqÞ@q

@x1

¼ 0 ð31Þ

with f1and B1deﬁned by(18)and the Riemann data QLand QR.

When an approximate Riemann solver is used for the numerical resolution, this would result in three propagating discon-
tinuities that are moving with speeds k^{1;k}_{i1=2;j}and the jumps W^{1;k}_{i1=2;j}across each of them for k = 1, 2, 3, see[31,33,34]for an
example.

(3) Deﬁne scaled speeds
k^{1;k}_{i1=2;j}¼ ðS2Þ_{i1=2;j}k^{1;k}_{i1=2;j}

and rotate jumps back to the Cartesian coordinates by
W^{1;k}_{i1=2;j}¼ R^{T}_{i1=2;j}W^{1;k}_{i1=2;j}

for k = 1, 2, 3.

(4) Determine the left- and right-moving ﬂuctuations in the form

A^{}_{1}DQ_{i1=2;j}¼X^{3}

k¼1

k^{1;k}_{i1=2;j}

W^{1;k}_{i1=2;j}:

As usual, the notations for the quantities k^{±}are set by k^{+}= max(k, 0) and k^{}= min(k, 0).

In a similar manner, we may determine the up- and down-moving ﬂuctuations at the edge (i, j 1/2) in the form

A^{}2DQ_{i;j1=2}¼X^{3}

k¼1

k^{2;k}_{i;j1=2}

W^{2;k}_{i;j1=2}

that is as the result of solving

@q

@tþ1 J

@f2ðqÞ

@n2

þ B2ðqÞ1 J

@q

@n2

¼ 0
with the initial data Q^{n}_{i;j1}and Q^{n}_{i;j}.

3.2. High-resolution corrections

To achieve high-resolution (i.e., second order accurate on smooth solutions, and sharp and monotone proﬁles on discon- tinuous solutions), the speeds and the limited version of the jumps are used to construct the piecewise linear correction terms as before (cf.[21]), and are added to(30)in ﬂux difference form as

Q^{nþ1}_{ij} :Q^{nþ1}_{ij} 1

### j

ijDt Dn1

Fe^{1}_{iþ1=2;j} eF^{1}_{i1=2;j}

1

### j

ijDt Dn2

Fe^{2}_{i;jþ1=2} eF^{2}_{i;j1=2}

: Here at the edge (i 1/2, j) the correction ﬂux takes the form

Fei1=2;j¼1 2

X^{3}

k¼1

k^{1;k}_{i1=2;j}

1 Dt

### j

i1=2;jDn1k^{1;k}_{i1=2;j}

Wf^{1;k}_{i1=2;j};

and analogously at the edge (i, j 1/2) the correction ﬂux has the form

Fe^{2}_{i;j1=2}¼1
2

X^{3}

k¼1

k^{2;k}_{i;j1=2}

1 Dt

### j

_{i;j1=2}Dn2

k^{2;k}_{i;j1=2}

Wf^{2;k}_{i;j1=2};

where

### j

i1/2,j= (### j

i1,j+### j

i,j)/2 and### j

i,j1/2= (### j

i,j1+### j

i,j)/2. The quantity fW^{m;k}is a limited value of W

^{m;k}obtained by comparing W

^{m;k}with the corresponding W

^{m;k}from the neighboring Riemann problem to the left (if k

^{m,k}> 0) or to the right (if k

^{m,k}< 0) for m = 1, 2 and k = 1, 2, 3.

In addition to that, a transverse propagation of wave is also included in the method as a part of the high-resolution cor-
rection terms. Here the right-moving ﬂuctuation A^{þ}_{1}DQ_{i1=2;j}, for instance, is decomposed into transverse ﬂuctuations
A^{}_{2}A^{þ}_{1}DQi1=2;jwhich can be used to update the solutions above and below cell (i, j). In a similar manner, the left-moving ﬂuc-
tuation A^{}_{1}DQi1=2;jis split into A^{}_{2}A^{}_{1}DQi1=2;jwhich are used to update the solutions above and below cell (i 1, j), see[21]

for the details.

With the transverse wave-propagation, the method is typically stable as long as the time stepDt satisﬁes a variant of the CFL (Courant–Friedrichs–Lewy) condition of the form

### m

¼Dt maxi;j;k

k^{1;k}_{i1=2;j}

### j

ip;jDn1;
k^{2;k}_{i;j1=2}

### j

i;jpDn20

@

1

A 6 1; ð32Þ

where ip= i if k^{1;k}_{i1=2;j}>0 and i 1 if k^{1;k}_{i1=2;j}<0; jpis deﬁned analogously (cf.[40]). Furthermore, by following the basic steps
discussed in[9,21], it can be shown that the method is quasi-conservative in the sense that when applying the method to our
multiphase model(15)not only the conservation laws but also the transport equations are approximated in a consistent
manner by the method.

3.3. Three-dimensional extension

To extend this mapped grid method from two to three space dimensions, we use hexahedral meshes in place of quadri- lateral-grid cells, seeFig. 11for an example. We use the ﬂuctuation form of the method as usual in three dimensions for the solution updates that is an easy generalization of the three-dimensional wave-propagation method on Cartesian grids pro- posed by Langseth and LeVeque[19]and also the two-dimensional method on mapped grid[21]. To implement the method, it is useful to make reference to the CLAWPACK webpage, see[9,36]in particular, for the programming details. In Section4.2, we present some sample results that show the feasibility of this method to practical compressible multiphase problems.

4. Numerical results

We now present numerical results to validate our mapped grid algorithm for compressible multiphase ﬂow problems in two and three dimensions. Note that in this section we have only present solutions obtained using the

### c

-based model(13)to the method. This is because we have found a little difference between the results as compared to the ones using the### a

-based model(27)to the method for simulations.4.1. Two-dimensional case

4.1.1. Smooth vortex ﬂow

We begin our tests by performing a convergence study of the computed solutions for a two-dimensional vortex evolution problem (cf.[12,29,43]) that shows the order of accuracy that is attained for our high-resolution method as the mesh is re- ﬁned. In this problem, we assume a single-phase ideal gas ﬂow with

### c

= 1.4 and B ¼ 0 in the stiffened gas equation of state (2). Initially, over a square domain of size [0, 10] [0, 10], the state variables for the vortex are set by### q

¼ 1 ð### c

1Þ^{2}

8

### cp

^{2}

^{expð1 r}

2Þ

1=ðc1Þ

; p ¼

### q

^{c};

u1¼ 1

2

### p

^{expðð1 r}

2Þ=2Þ xð 2 x2Þ;

u2¼ 1 þ

2

### p

^{expðð1 r}

2Þ=2Þðx1 x1Þ;

where r ¼ ﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
ðx1 x_{1}Þ^{2}þ ðx2 x_{2}Þ^{2}
q

is the distance between points (x_{1}, x_{2}) and the vortex center ðx_{1}; x_{2}Þ ¼ ð5; 5Þ, and

### q

, u1, and u2to the mean ﬂow with### q

¼ 1; p ¼ 1; u1¼ 1, and u2¼ 1, and it is known that with the periodic conditions in both the x1- and x2-direction the exact solution for this problem is simply the passive motion of a smooth vortex by the mean ﬂow velocity.We compare the behavior of the high-resolution wave-propagation method described in Section3on three different types of grids as illustrated inFig. 2:

Grid 1: Cartesian grids with square grid cells.

Grid 2: Quadrilateral grids of the type described in[10].

Grid 3: Quadrilateral grids of the type described in[9].

see[22]for a similar computation on the transport of a scalar quantity in a divergence-free velocity ﬁeld. In order to estimate the order of accuracy of the method, for each grid type, we compute the solution at the mesh reﬁnement sequence:

Dn^{ðiÞ}_{1} ¼Dn^{ðiÞ}_{2} ¼Dn^{ðiÞ}¼ 1=2^{ðiþ2Þ}for i = 0, 1, 2, 3, i.e., by using an N N grid for N = 40, 80, 160 and 320. To ensure that our conver-
gence study is not adversely affected by a loss of accuracy near local extrema in the solution, there is no limiter being used in
the computations. In addition, the Courant number

### m

= 0.9 deﬁned by(32), and the Roe approximate Riemann solver are employed in the tests.Numerical results at time t = 10 obtained from all the 12 tests are presented inTables 1–3. Here E1ðzÞ and EmðzÞ denote in turn the sequence of the one- and maximum-norm error of the cell-average solution in z to the true solution at the cell center for z =

### q

, u_{1}, u

_{2}, and p. We estimate the rate of convergence using the errors on two consecutive grids based on the formula

convergence order ¼ ln E^{ði1Þ}_{k} ðzÞ
E^{ðiÞ}_{k}ðzÞ
!,

ln Dn^{ði1Þ}
Dn^{ðiÞ}

!

for k = 1, m and i = 1, 2, 3.

From the tables, we observe that the method is second order accurate in most cases on this problem. It is important to note that the error behavior of the method on Grid 2 is on the same order of magnitude on Grid 1, the Cartesian grid, while this is not the case on Grid 3 which is a less smooth grid as compared to Grid 2. This indicates that the use of a smooth grid in a mapped grid method can indeed give a better resolution of the solution.

0 5 10

0 2 4 6 8 10

0 5 10

0 2 4 6 8 10

0 5 10

0 2 4 6 8 10

Fig. 2. Three types of grids used for the smooth vortex ﬂow problem.

Table 1

High-resolution results for the smooth vortex test on Grid 1; one- and maximum-norm errors in primitive variables are shown.

N E1ðqÞ Order E1ðu1Þ Order E1ðu2Þ Order E1ðpÞ Order

40 0.6673 2.3443 1.7121 0.8143

80 0.1792 1.90 0.6194 1.92 0.4378 1.97 0.2128 1.94

160 0.0451 1.99 0.1537 2.01 0.1104 1.99 0.0536 1.99

320 0.0113 2.00 0.0384 2.00 0.0276 2.00 0.0134 2.00

N EmðqÞ Order Emðu1Þ Order Emðu2Þ Order EmðpÞ Order

40 0.1373 0.3929 0.1810 0.1742

80 0.0377 1.87 0.1014 1.95 0.0502 1.85 0.0482 1.85

160 0.0093 2.02 0.0248 2.03 0.0123 2.03 0.0119 2.02

320 0.0022 2.07 0.0062 2.00 0.0030 2.04 0.0029 2.04

Table 2

High-resolution results for the smooth vortex test on Grid 2; one- and maximum-norm errors in primitive variables are shown.

N E1ðqÞ Order E1ðu1Þ Order E1ðu2Þ Order E1ðpÞ Order

40 0.9298 2.6248 2.1119 1.2104

80 0.2643 1.81 0.7258 1.85 0.5296 2.00 0.3277 1.89

160 0.0674 1.97 0.1833 1.99 0.1309 2.02 0.0845 1.96

320 0.0169 2.00 0.0458 2.00 0.0327 2.00 0.0212 1.99

N EmðqÞ Order Emðu1Þ Order Emðu2Þ Order EmðpÞ Order

40 0.1676 0.4112 0.2259 0.2111

80 0.0471 1.83 0.1242 1.73 0.0645 1.79 0.0586 1.85

160 0.0126 1.91 0.0333 1.90 0.0162 2.02 0.0149 1.97

320 0.0033 1.93 0.0085 1.97 0.0040 2.00 0.0038 1.98

4.1.2. Passive interface evolution

We are next concerned with an interface-only problem that the exact solution consists of a circular water column evolv-
ing in the air with uniform equilibrium pressure p ¼ 10^{5}Pa and constant particle velocity (u_{1}; u_{2}Þ = (10^{3}, 10^{3}) m/s throughout
a quarter annulus domain. Here we take the initial condition that inside the column of radius r0= 0.2 m about the center
ðx1; x2Þ ¼ ð0:8; 0:8Þ m the ﬂuid is water with the data

ð

### q

;### c

;### q

0;BÞr6r_{0}¼ ð10

^{3}kg=m

^{3};4:4; 10

^{3}kg=m

^{3};2:64 10

^{6}ðm=sÞ

^{2}Þ;

while outside the column the ﬂuid is air with the data

ð

### q

;### c

;### q

_{0};BÞr>r

_{0}¼ ð1:2 kg=m

^{3};1:4; 1:2; 0Þ:

Note that despite the simplicity of the solution structure, this problem is one of the popular tests for the numerical validation of a compressible multiphase ﬂow solver (cf.[30,43]).

Table 3

High-resolution results for the smooth vortex test on Grid 3; one- and maximum-norm errors in primitive variables are shown.

N E1ðqÞ Order E1ðu1Þ Order E1ðu2Þ Order E1ðpÞ Order

40 4.8272 4.7734 5.3367 5.4717

80 1.5740 1.62 1.5633 1.61 1.5660 1.77 1.5634 1.81

160 0.4536 1.79 0.4559 1.78 0.4537 1.79 0.4560 1.78

320 0.1215 1.90 0.1221 1.90 0.1222 1.89 0.1221 1.90

N EmðqÞ Order Emðu1Þ Order Emðu2Þ Order EmðpÞ Order

40 0.4481 0.4475 0.4765 0.4817

80 0.1170 1.94 0.1181 1.92 0.1196 1.99 0.1191 2.02

160 0.0434 1.43 0.0431 1.45 0.0442 1.43 0.0440 1.44

320 0.0117 1.89 0.0119 1.86 0.0119 1.89 0.0118 1.89

Fig. 3. Numerical results for an interface-only problem in a quarter annulus at time t = 520ls. On the top, surface plots of the density and pressure are shown, and on the bottom, cross-sectional plot (along the line n2=p/4) of the density and scatter plots of the relative error of the pressure are displayed, respectively. Here the solid line is the exact solution, and the dashed line is the initial condition at the time t = 0.

To discretize this quarter-annulus region, we use polar coordinates:

x1ðn1;n2Þ ¼ n1cosðn2Þ;

x2ðn1;n2Þ ¼ n1sinðn2Þ;

for 0.5 m 6 n162.5 m and 0 6 n26

### p

/2, see Chapter 23 of[21]for an illustration.Fig. 3shows numerical results of the den- sity and pressure at time t = 520### l

s obtained using our algorithm with theMINMODlimiter and the Roe solver on a 100 100 polar grid. From the 3D surface plot and the cross-section plot (along n2=### p

/4) of the density, we observe that the water col- umn retains its circular shape and appears to be very well located also. From the 3D surface plot of the pressure and the scatter plot the relative error of the pressure, we ﬁnd the computed pressure remains in the correct equilibrium state p (to be more accurate, the difference of these two is only on the order of machine epsilon), without any spurious oscillations near the air–water interface. Here, we use non-reﬂecting boundary conditions on all sides of the quarter annulus, while car- rying out the computations.4.1.3. Moving cylindrical vessel

Our next example concerns a moving cylindrical vessel problem studied by Banks et al.[5]in that inside the circle of ra- dius r0= 0.8 about the origin, there is a planar material interface located initially at x1= 0 that separates air on the left with the state variables

ð

### q

;u1;u2;p;### c

;### q

0;BÞ ¼ ð1; 1; 0; 1; 1:4; 1; 0Þ;and helium on the right with the state variables

ð

### q

;u1;u2;p;### c

;### q

_{0};BÞ ¼ ð0:138; 1; 0; 1; 1:67; 0:138; 0Þ:

Note that in this set up we have imposed a uniform ﬂow velocity (u1, u2) = (1, 0) throughout the domain, and so we are in the frame of the vessel moving with speed one in the x1-direction.

To ﬁnd an approximate solution of this problem, we use a mapped-grid approach proposed by Calhoun et al.[9]in that a grid point (n1, n2) in the computational domain [1, 1] [1, 1] is mapped to a grid point (x1, x2) in the circular domain by some simple algebraic rules, seeFig. 4for an illustration. Numerical Schlieren images and pseudo-color plots of pressure are shown inFig. 4at four different times t = 0.25, 0.5, 0.75, and 1.0, where a 800 800 grid is used in the run. From the ﬁg- ure, it is easy to see that due to the impulsive motion of the vessel a rightward-going shock wave and a leftward-going rar- efaction wave emerge from the left- and right-side boundary, respectively. Subsequently, these two waves would be interacting with the material interface that leads to collision of various transmitted and reﬂected waves. When we compare our results with those ones appeared in the literature (cf.[5,43]), as far as the global wave structures are concerned, we no- tice good qualitative agreement of the solutions.

To check the quantitative information of our computed solutions,Fig. 5compares the cross-sectional results for the same run along the circular boundary with those obtained using a wave-propagation based Cartesian grid embedded boundary method (cf.[20,35,37]). The close agreement between the mapped grid and Cartesian grid results in the shock waves and material interfaces are clearly observed.

4.1.4. Shock–bubble interaction in a nozzle

As an example to show how our algorithm works on shock waves in a more general two-dimensional geometry, we are interested in a shock–bubble interaction problem in a nozzle. For this problem, the shape of the nozzle is described by a ﬂat curve

x^{t}_{2}¼ 1 m

on the top, and by the witch of Agnesi

x^{b}_{2}ðx1Þ ¼ 8a^{3}
x^{2}_{1}þ 4a^{2}

on the bottom for a = 0.2 m and 2 m 6 x_{1}63 m. We use the initial condition that is composed of a planarly rightward-
going Mach 1.422 shock wave located at x1= 1.8 m in liquid traveling from left to right, and a stationary gas bubble of
radius r0= 0.2 m and center ðx1; x2Þ ¼ ð1; 0:5Þ m in the front of the shock wave. Inside the gas bubble, we have the data

ð

### q

;u1;u2;pÞ ¼ ð1:2 kg=m^{3};0; 0; 10

^{5}PaÞ;

while outside the gas bubble where the ﬂuid is liquid, we have the preshock state ð

### q

;u1;u2;pÞ ¼ ð10^{3}kg=m

^{3};0; 0; 10

^{5}PaÞ;

and the postshock state

ð

### q

;u1;u2;pÞ ¼ ð1:23 10^{3}kg=m

^{3};432:69 m=s;0; 10

^{9}PaÞ:

Here the material-dependent parameters ð

### c

;### q

_{0};BÞ for the gas- and liquid-phase are taken as (1.4,1.2 kg/m

^{3}, 0) and (4.4, 10

^{3}kg/m

^{3}, 2.64 10

^{6}(m/s)

^{2}), respectively.

Fig. 4. Numerical Schlieren images (on the left) and pseudo colors of pressure (on the right) for an impulsively driven cylinder containing an air–helium material interface. Solutions from top to bottom are at times t = 0.25, 0.5, 0.75, and 1.0.

In carrying out the computation, we consider a body-ﬁtted quadrilateral grid with the mapping function

x1ðn1;n2Þ ¼ n1;

x2ðn1;n2Þ ¼ x^{b}_{2}ðn1Þ n^{t}_{2} n2

n^{t}_{2} n^{b}_{2}

!

for 2 m 6 n163 m, 0 6 n2 61 m;n^{b}_{2}¼ 0, and n^{t}_{2}¼ 1 m. The boundary conditions are the supersonic inﬂow on the left-hand
side, the non-reﬂecting on the right-hand side, and the solid wall on the remaining sides. InFig. 6, we show the Schlieren
images and pseudo colors of pressure at six different times t ¼ 0:3; 0:5; 0:7; 1:2; 1:6, and 2.5 ms obtained using a
1000 200 grid. From the ﬁgure, it is easy to see that after the passage of the shock to the gas bubble, the upstream wall
begins to spall across the bubble, yielding a refracted air shock traveling within it until its ﬁrst reﬂection on the downstream
bubble wall. Noticing that this upstream bubble wall would be involute eventually to form a jet which subsequently crosses
the bubble and sends an intense blast wave out into the surrounding liquid. In the meantime, the incident shock wave along
the bottom curved boundary would be diffracted into a simple Mach reﬂection.

To see the quantitative information of the solutions,Fig. 7compares the cross-section of the results along the n_{2}= 0.5 m
line with the same method but with a ﬁner 2000 400 grid. Convergence of our computed solutions to the correct weak
ones is clearly observed. In addition,Fig. 8 compares the cross-section of the results for the same run along the bottom
boundary with the results obtained using a Cartesian grid embedded boundary method[35]. The close agreement between
the mapped grid and Cartesian grid results are seen also.

4.1.5. Underwater explosion with circular obstacles

Our next example for problems in more general geometries is an underwater explosion ﬂow with circular obstacles, see
[24]for a similar computation but with a square obstacle. In this test, the physical domain is a rectangular region of size
([2, 2] [1.5, 1]) m^{2}in which inside the domain there are two circular obstacles, denoted by D1and D2, with the centers
ðx^{D}_{1}^{1}; x^{D}_{2}^{1}Þ ¼ ð0:6; 0:8Þ m and ðx^{D}_{1}^{2}; x^{D}_{2}^{2}Þ ¼ ð0:6; 0:4Þ m, respectively, and of the same radius rD1¼ rD2¼ 0:2 m. The initial
condition we consider is composed of a horizontal air–water interface at x2= 0 and a circular gas bubble in water that lies
below the interface. Here all the ﬂuid components are at rest initially. When x2P0, the ﬂuid is air with the state variables
Fig. 5. Cross-sectional plots of the results for the moving vessel run shown inFig. 4along the circular boundary, where the solid lines are results obtained
using a Cartesian grid embedded boundary method[35]with the same grid size.