Journal of Computational Physics

(1)

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 efficient numerical simulation of compressible multiphase flow 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 fluid mixing when approximating 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 flow problems. We validate our algorithm by performing numerical tests in two and three dimensions that show second order accurate results for smooth flow 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.

1. Introduction

Our goal is to describe a simple mapped-grid approach for efficient numerical resolution of compressible multiphase flow in general multi-dimensional geometries. As a first endeavor towards the method development, we are concerned with a simplified model problem, where the flow 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 fluid 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 fluid component in a Cartesian coordinates can be written as

@

@t

q q

ui

E 0 B@

1 CA þX^N^d

j¼1

@

@xj

q

uj

q

uiujþ pdij

Eujþ 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.

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

(2)

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

₀Þ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 þPNd

j¼1

q

u²_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 flows on mapped grids to the case of a multiphase flow. It is well known that the principal problem in the usual extension is the occurrence of spurious pressure oscillations when two or more fluid components are present in a grid cell (cf.[5]and references therein). Here the algorithm uses a curvilinear coordinate formulation of a fluid-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-dependent material quantities. In this approach, as in its Cartesian coordinates counterpart (cf.[30,31]), the latter equations are derived to ensure the correct fluid 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 significant 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-fitted mapped grid than on a fixed 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 efficient 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 simplified homogeneous multiphase flow in curvilinear coordinates. In Section3, we review briefly 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 fundamental 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 J

X^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Þ

(3)

where Uj¼PN_d

i¼1uiJ_jiis 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₁₁ J₁₂ J₁₃ J₂₁ J₂₂ J₂₃ J₃₁ J₃₂ J₃₃ 0

B@

1 CA ¼

b2c3 b3c2 a3c2 a2c3 a2b3 a3b2

0 B@

1

CA ð5Þ

and the quantity J = detj@(x₁, x₂, x₃)/@(n₁, n₂, n₃)j is the Jacobian of the mapping which can be computed by

J ¼X³

i¼1

aiJ_1i¼X³

i¼1

biJ_2i¼X³

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 equations

To 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 constant 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 J

X^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

₀

c

1 B

þ1 J

X^N^d

j¼1

Uj

@

@nj

0, and B.

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

₀B

c

1

þ1 J

X^N^d

j¼1

1

þ1 J

X^N^d

j¼1

@

@nj

q

B

c

1Uj

¼ 0; ð12Þ

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

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Þ²

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

1

T

;

; 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

¼ 0; ð17Þ

with fjand Bjdeﬁned in turn by

; 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.

(5)

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_jiis 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 equations

_0;i

c

_i 1 Bi

¼ p

c

1

q

₀

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

p

c

1¼X^M^f

¼P^Mf

i¼1

a

i

q

_i, we are able to rewrite them straightforwardly into a componentwise form as

@

@t

a

i

c

_i 1

_0;iBi

c

i 1

þ1 J

X^N^d

j¼1

Uj

@

@nj

a

i

q

_0;iBi

c

i 1

¼ 0; ð24Þ

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

c

q

_iÞ þ1 J

X^N^d

j¼1

@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 equation of state

p ¼ E PN_d

i¼1ð

q

uiÞ²

2

q

^þ

_i 1; where we have assumed

a

Mf ¼ 1 PMf1

i¼1

a

i.

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

To end this section, we comment that if x₁is the axisymmetric direction, an axisymmetric version of our multiphase models in two dimensions can be written as

@q

@tþ1 J

X²

j¼1

@

@nj

¼ wðqÞ; ð28Þ

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

; ð29Þ

when the

c

-based model is considered, and have

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

T

;

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-

(7)

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ⁿ_ij 1 MðCijÞ

Z

Cij

qðx1;x2;tnÞdx1dx2¼ 1

j

ijDn1Dn2

Z bCij

ij

Dt Dn2

A^þ2DQ_i;j1=2þ A2DQ_i;jþ1=2

: ð30Þ

Here A^þ₁DQi1=2;j;A₁DQiþ1=2;j;A^þ₂DQi;j1=2, and A₂DQi;jþ1=2 are the right-, left-, up-, and down-moving ﬂuctuations, respectively, 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₁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ⁿ_i1;j; if n1<ðn1Þ_i1=2; Qⁿ_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;â2Þ_i1=2;jand~ti1=2;j¼ ðâ2; ^b2Þ_i1=2;jbe the unit normal and tangential vectors to the cell edge (i 1/2, j) in the physical grid, where âi¼ ai=Si and ^bi¼ bi=Si are the scaled version of the metric elements ai and bi in (7) with Si¼

ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi a²_i þ b²_i q

for i = 1, 2. Then to compute A₁DQ_i1=2;j, as in[9,21], we perform the following steps:

(1) Transform the data Qⁿ_i1;jand Qⁿ_i;jinto the new data QLand QRvia QL¼ Ri1=2;jQⁿ_i1;j; QR¼ Ri1=2;jQⁿ_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).

(8)

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 discontinuities that are moving with speeds k^1;k_i1=2;jand the jumps W^1;k_i1=2;jacross 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;jk^1;k_i1=2;j

and rotate jumps back to the Cartesian coordinates by W^1;k_i1=2;j¼ R^T_i1=2;jW^1;k_i1=2;j

for k = 1, 2, 3.

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

A₁DQ_i1=2;j¼X³

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

A2DQ_i;j1=2¼X³

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

_i;j1=2Dn2

k^2;k_i;j1=2

Wf^2;k_i;j1=2;

i,j)/2. The quantity fW^m;kis a limited value of W^m;kobtained by comparing W^m;kwith the corresponding W^m;kfrom 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.

(9)

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 fluctuation A^þ₁DQ_i1=2;j, for instance, is decomposed into transverse fluctuations A₂A^þ₁DQi1=2;jwhich can be used to update the solutions above and below cell (i, j). In a similar manner, the left-moving fluctuation A₁DQi1=2;jis split into A₂A₁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 max

i;j;k

k^1;k_i1=2;j

j

ip;jDn1

; k^2;k_i;j1=2

j

i;jpDn2

0

@

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 quadrilateral-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 proposed 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Þ

²

8

cp

² ^{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 ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðx1 x₁Þ²þ ðx2 x₂Þ² q

is the distance between points (x₁, x₂) and the vortex center ðx₁; x₂Þ ¼ ð5; 5Þ, and

= 5 is the vortex strength. Note that the above states are as the results of an isentropic perturbation in p/

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:

(10)

Dn^ðiÞ₁ ¼Dn^ðiÞ₂ ¼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 convergence 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₁, u₂, 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

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

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

(11)

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⁵Pa and constant particle velocity (u₁; u₂Þ = (10³, 10³) 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

Table 3

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

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.

(12)

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 density 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 column 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 carrying 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 radius 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

To find 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 figure, 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 reflected 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₂¼ 1 m

on the top, and by the witch of Agnesi

x^b₂ðx1Þ ¼ 8a³ x²₁þ 4a²

on the bottom for a = 0.2 m and 2 m 6 x₁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

₀;BÞ for the gas- and liquid-phase are taken as (1.4,1.2 kg/m³, 0) and (4.4, 10³kg/m³, 2.64 10⁶(m/s)²), 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.

(14)

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₂ðn1Þ n^t₂ n2

n^t₂ n^b₂

!

for 2 m 6 n163 m, 0 6 n2 61 m;n^b₂¼ 0, and n^t₂¼ 1 m. The boundary conditions are the supersonic inflow on the left-hand side, the non-reflecting 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 figure, 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 first reflection 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 reflection.

To see the quantitative information of the solutions,Fig. 7compares the cross-section of the results along the n₂= 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 flow 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²in which inside the domain there are two circular obstacles, denoted by D1and D2, with the centers ðx^D₁¹; x^D₂¹Þ ¼ ð0:6; 0:8Þ m and ðx^D₁²; x^D₂²Þ ¼ ð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 fluid components are at rest initially. When x2P0, the fluid 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.