## A wave-propagation based volume tracking method

## for compressible multicomponent ﬂow in two space dimensions

### Keh-Ming Shyue

^{*}

Department of Mathematics, National Taiwan University, Taipei 106, Taiwan Received 18 June 2005; received in revised form 25 October 2005; accepted 27 October 2005

Abstract

We present a simple volume-of-ﬂuid approach to interface tracking for inviscid compressible multicomponent ﬂow problems in two space dimensions. The algorithm uses a uniform Cartesian grid with some grid cells subdivided by tracked interfaces, approximately aligned with the material interfaces in the ﬂow ﬁeld. A standard volume-moving procedure that consists of two basic steps: (1) the update of a discrete set of volume fractions from the current time to the next and (2) the reconstruction of interfaces from the resulting volume fractions, is employed to ﬁnd the new location of the tracked interfaces in piecewise linear form at the end of a time step. As in the previous work by LeVeque and the author to front tracking based on a surface-moving procedure (R.J. LeVeque, K.-M. Shyue, Two-dimensional front tracking based on high-resolution wave propagation methods, J. Comput. Phys. 123 (1996) 354–368), a conservative high-resolution wave propagation method is applied on the resulting slightly non-uniform grid to update all the cell values, while the stability of the method is maintained by using a large time step idea even in the presence of small cells and the use of a time step with respect to the uniform grid cells. We validate our algorithm by performing the simulation of a Mach 1.22 shock wave in air over a circular R22 gas bubble, where sensible agreement of some key ﬂow features of the computed solutions are observed when direct comparison of our results are made with the existing experimental and numerical ones that appear in the lit- erature. Other computations are also presented that show the feasibility of the algorithm together with a mixture type of the model equations developed by the author (K.-M. Shyue, A ﬂuid-mixture type algorithm for compressible multicom- ponent ﬂow with Mie–Gru¨neisen equation of state, J. Comput. Phys. 171 (2001) 678–707) for practical multicomponent problems with general compressible materials characterized by a Mie–Gru¨neisen form of the equation of state.

2005 Elsevier Inc. All rights reserved.

AMS: 65M06; 76L05; 76M20; 76T05

Keywords: Volume tracking; Wave propagation method; Multicomponent ﬂows; Mie–Gru¨neisen equation of state; Impact problem

0021-9991/$ - see front matter 2005 Elsevier Inc. All rights reserved.

doi:10.1016/j.jcp.2005.10.030

* Tel.: +8862 3366 2866; fax: +8862 2391 4439.

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

Journal of Computational Physics xxx (2005) xxx–xxx

www.elsevier.com/locate/jcp

1. Introduction

Our goal in this paper is to describe a simple volume tracking approach for the eﬃcient numerical resolu- tion of material interfaces arising from inviscid compressible multicomponent problems in two space dimen- sions. We use a Eulerian formulation of the equations in which the principal motion of the single-component ﬂow is governed by the basic conservation laws for the mass, momenta in the x- and y-direction, and energy as

o ot

q qu qv qE 0 BB B@

1 CC CAþ o

ox

qu
qu^{2}þ p

quv qEuþ pu 0

BB B@

1 CC CAþ o

oy

qv
quv
qv^{2}þ p
qEvþ pv
0

BB B@

1 CC

CA¼ 0; ð1Þ

respectively. Here q denotes the density, u is the particle velocity in the x-direction, v is the particle velocity in the y-direction, p is the pressure, and E is the speciﬁc total energy. To complete the model, the constitutive law for the thermodynamic behavior of the ﬂuid component of interests is assumed to satisfy a Mie–Gru¨neisen equation of state of the form

pðq; eÞ ¼ p_{ref}ðqÞ þ CðqÞ qe qe½ refðqÞ; ð2Þ

where e represents the speciﬁc internal energy, C = (1/q)(op/oe)jqis the Gru¨neisen coeﬃcient, and pref, erefare
the properly chosen states of the pressure and internal energy along some reference curve (e.g., along an
isentrope, a single shock Hugoniot, or other empirically ﬁt curves) in order to match the experimental data
of the material being examined, see [54] and references cited therein for more details. As usual, we have
E = e + (u^{2}+ v^{2})/2.

As in the previous work by LeVeque and the author[31]to general front tracking for nonlinear hyperbolic systems of conservation laws in two dimensions, our approach to volume tracking is based on a fully discret- ized version of the high-resolution ﬁnite volume method in wave-propagation form. The algorithm uses an underlying uniform Cartesian grid with some rectangles subdivided by the tracked interfaces into two parts, where discontinuities in the material interfaces (or slip lines for that matter) are expected. Rather than employ- ing a marker-and-cell representation of the interface and a surface-moving procedure to advance the interface location as before, which can be done for some simple cases (cf.[51,52]) but is very troublesome (at least in view of the code development aspect) in dealing with situations such as the splitting or merging of interfaces and the changes of complex multidimensional interfacial topology (cf.[9]), we adopt a volume-fraction based interface-moving procedure for an easy treatment of the aforementioned issues instead, see [19,38] for the other possible alternates. That is to say, we introduce a scalar ﬁeld Y in the algorithm that is deﬁned to be the fraction of each cell that is on one particular side of the interface. In the current interest of two-dimen- sional compressible ﬂow problems, this volume-fraction function Y would be updated in each time step by solving an evolution equation of the form

oY ot þ uoY

oxþ voY

oy ¼ 0 ð3Þ

with a standard wave-propagation method. Here u and v are the particle velocities governed by the full Euler equations(1) which is unlike the case for many incompressible ﬂow problems where a divergence-free con- straint should be fulﬁlled on them for the fundamental agreement of the mathematical model being used.

When this is done, the new interface can then be reconstructed in a cell-by-cell fashion from the information on Y, yielding a set of disconnected line segments (this is as opposed to a connected line segment case obtained by using a typical surface-moving procedure[6,31]) in the cells containing parts of the interface, see Section3 for more exposition.

An important property of the wave-propagation form of the method is that reasonable time steps can be taken even if some of the subcells created by the tracked interfaces are orders of magnitude smaller than the underlying Cartesian cells. Uniform time steps are used throughout the computation, with the time step cho- sen so that the usual CFL (Courant–Friedrichs–Levy) condition is satisﬁed relative to the size of the regular grid cells. The method would typically remain stable in this case because near the small subcells the numerical

domain of dependence is extended automatically based on a large time step technique developed by LeVeque for general hyperbolic systems of partial diﬀerential equations (cf. [24–27]). High-resolution results can be achieved in a straightforward manner in most part of the domain, except possibly for cells near the tracked interfaces. This is discussed further in Section2.

It should be mentioned that as in our surface-based front tracking algorithm[31], the basic philosophy of the proposed volume tracking method is that the subdivision of cells is not assumed to give the deﬁnitive loca- tion of the true interface, but is viewed rather as an approximate location yielding a reﬁned grid that is better able to represent the discontinuous solution than the Cartesian grid alone. As a result, there is minimal smear- ing of the interface in the solution. To prevent the occurrence of spurious oscillation in the pressure when such a scenario occurs in the multicomponent context, the algorithm uses a mixture type of the model system described in[54]as the basis, and a standard ﬁnite volume method based on the wave-propagation formula- tion as for the consistent and accurate approximation of the system. Numerical results presented in Section5 show that this gives an eﬃcient way for the improvement of the resolution of the interfaces, and the correct convergence behavior of the weak solution for shock waves.

It is without doubt that there are many instances to have the desire of keeping the interface sharp, and con- sequently the advantage of a tracking method over a ﬁxed grid capturing method. One good example among them is in the simulation of compressible two-phase ﬂows with very diﬀerent ﬂuid component separated by interfaces, and under extreme ﬂow condition, see[18,32,37,65,66,69]for an example. In this case, despite some simple situations (cf.[54]), numerical experiences indicate that it would be very diﬃcult in general to devise a pressure-oscillation-free interface-capturing method to deal with the occurrence of more than one ﬂuid com- ponent within a grid cell. In fact, it is one of the main motivations of the present attempt to develop a robust volume tracking method that is viable for a wide variety of two-phase ﬂow problems with real materials.

There have been quite a few eﬀorts over the years to develop eﬃcient volume-of-ﬂuid interface tracking methods for various practical applications. In regards to incompressible ﬂow problems governed by the Navier–Stokes equations, for example, representative works in this direction are such as Hirt and Nichols for free boundary problems [16], Welch and Wilson for phase transition problems[67], and Gueyﬃer et al.

[13], Lafaurie et al. [22], Rudman [47,48], and Sussman and Puckett [59] for two-phase ﬂow problems. In regards to the current compressible ﬂow problems governed by the inviscid Euler equations, there are works done by Chern and Colella[5]for nonlinear hyperbolic systems of conservation laws (see[1]also for a version of the method with the local adaptive mesh reﬁnement included), and Miller and Puckett[37]for problems in multiple condensed matters.

The approach we take here is in essence similar to that of Chern and Colella [5]in which they also use a volume-of-ﬂuid approach to keep track of the interface location, and then apply high-resolution conser- vative ﬁnite volume methods on the resulting grid. The key diﬀerence is that they use a more traditional ﬂux diﬀerence form, together with a ﬂux redistribution algorithm to maintain stability in the presence of small cells. It has been demonstrated before (cf. [31,52]) that our wave-propagation form can give good results when it is combined with a surface-based front tracking algorithm in one and two dimensions.

One major objective of this paper is to show that this method can still be very eﬀective when it is used in conjunction with a volume-based interface tracking method for two-dimensional compressible hydrody- namics problems.

The format of this paper is outlined as follows. In Section2, we review the basic idea of a ﬁnite volume method in wave-propagation form on a slightly non-uniform Cartesian grid with a moving interface. In Sec- tion3, we describe the standard piecewise linear interface construction method based on a given set of discrete volume fractions. The volume tracking algorithm for inviscid compressible ﬂow is discussed in Section 4.

Results of some sample validation tests as well as application of the method to practical problems are shown in Section5.

2. Wave-propagation methods

To set the groundwork for the later development of a volume-based interface tracking algorithm, we begin our discussion by reviewing a ﬁnite volume method in wave-propagation form for the numerical approxima- tion of general hyperbolic systems of conservation laws in two space dimensions,

oq

otþofðqÞ

ox þogðqÞ

oy ¼ 0. ð4Þ

Here q*2 R*^{m}denotes the vector of m conservative quantities, and f(q), g(q) denote the ﬂux vectors, see(1)for
an example. For all practical purposes as concerned below, it is essential to describe the method on a slightly
non-uniform grid that is composed of a set of regular cells with a ﬁxed mesh spacing Dx and Dy in the x- and
y-direction, respectively, and some irregular cells which are formed by inserting a piecewise linear representa-
tion of the tracked interfaces into the grid, yielding the subdivision of some regular cells into two subcells, see
Fig. 1a for an illustration.

In a ﬁnite volume formulation of the method, the cell average of each grid cell is deﬁned by integrating the solution over the cell and dividing by the area of the cell. Hence the approximate value of the cell average of the solution over any given cell S at time tncan be written as

Q^{n}_{S} 1
MðSÞ

Z

S

qðx; y; tnÞ dx dy; ð5Þ

where MðSÞ is the measure (area) of this cell. The methods we use are based on solving one-dimensional Rie- mann problems at each cell edge, and the waves (i.e., discontinuities moving at constant speeds) arising from the Riemann problem are employed to update the cell averages in the cells neighboring each edge (cf.

[26,27,31]).

2.1. First-order method

As a ﬁrst example to explain the basic idea of our ﬁrst-order method, we consider the cell edge between irregular cells C and D as illustrated inFig. 1a. We solve the one-dimensional Riemann problem normal to this face, which in this case will be

Fig. 1. First-order wave-propagation method. (a) The typical grid created by our volume tracking algorithm. (b) A wave propagating normally from the cell edge between cells C and D updates the solution in both cells D and F. (c) A wave propagating normally from the cell edge between cells C and E updates four neighboring cell averages. (d) A wave propagating both normally and tangentially from the cell edge between cells C and D updates the solution in cells D, F, and H. (e) A wave propagating both normally and tangentially from the cell edge between cells C and E updates ﬁve neighboring cell averages. Here Wprepresents the region swept out by the p-wave, and Wpq

represents the region swept out by the p-wave and q-wave in the normal and tangential to the cell edge, respectively, over a time step Dt for p,q = 1, 2, . . . , m.

oq

otþofðqÞ ox ¼ 0;

with the initial data given by Q^{n}_{C} and Q^{n}_{D}. If RoeÕs approximate solver is chosen to solve the Riemann problem,
the above nonlinear equation would be replaced by a linear system of the form

oq

otþ A Q ^{n}_{C}; Q^{n}_{D} oq
ox¼ 0;

where AðQ^{n}_{C}; Q^{n}_{D}Þ is a constant matrix that depends on the initial data and is a local linearization of the Jaco-
bian matrix of f(q) about an average state (cf.[29,46,62]). The result is that the jump Q^{n}_{D} Q^{n}_{C} is decomposed
into eigenvectors of the matrix A,

Q^{n}_{D} Q^{n}_{C}¼X^{m}

p¼1

apr_{p};

where rpis the pth eigenvector of A and Arp= kprpwith kpthe corresponding eigenvalue. The scalar apgives the wave strength.

Clearly, if the wave speed kp> 0 the pth wave propagates towards the right, while if kp< 0 it propagates towards the left.Fig. 1b shows a typical p-wave with kp> 0, where part of cell D and part of cell F are aﬀected by this wave over a time step Dt from the current time tnto the next tn+1. In the simplest case of the wave- propagation method, i.e., a two-dimensional extension of the ﬁrst-order GodunovÕs method, these cell aver- ages are updated by

Q^{nþ1}_{S} :¼ Q^{nþ1}_{S} MðWp\ SÞ

MðSÞ apr_{p} ð6Þ

for S = D, F. Here Wprepresents the region swept out by the wave and aprpis the jump across the wave. Note
that we have employed a standard initialization procedure that assigns each cell average to its values at the
previous time step, Q^{nþ1}_{S} :¼ Q^{n}_{S} for all S.

Now, if we are concerned with the cell edge between cells C and E (this is one segment of the tracked inter-
face), a Riemann problem would be solved in the direction normal to the interface with data Q^{n}_{C} and Q^{n}_{E}. Note
that the Euler equations (1)are isotropic and easily rotated to this frame (cf. [26]). The p-wave indicated in
Fig. 1c overlaps four cells and using (6) in a similar manner gives modiﬁcation of the cell average in each
of these cells by the jump across this wave, weighted by the fraction of the cell area covered by the wave.

To improve the stability of this ﬁrst-order method, it is a customary approach by introducing transverse propagation of waves based on the solutions of Riemann problems in the orthogonal direction (cf.

[26,27,31]). Again we use the cell edge between cells C and D inFig. 1a, as an example. Then for the equation we take the tangential portion of(4), i.e., the portion of equation in the y-direction,

oq

otþogðqÞ oy ¼ 0;

and for the initial data if we are working with the p-wave we take
Q^{L}¼ Q^{n}_{C}þX^{p1}

q¼1

aqr_{q}; Q^{R}¼ Q^{L}þ apr_{p}.

With the Roe solver, the solution to this Riemann problem yields a splitting of Q^{R} Q^{L}= aprpinto eigenvec-
tors of the Roe matrix B(Q^{L}, Q^{R}) of the form

apr_{p}¼X^{m}

q¼1

b_{q}w_{q}

where B is a suitable approximation of the Jacobian matrix of g(q), and B wq= lqwq. The wave shown in Fig. 1b is thus split into m waves propagating in the y-direction with speed given by the corresponding

eigenvalue lq, and a typical one is shown in Fig. 1d for the case lq> 0. Here to update the cell averages af- fected by this wave, we employ an analogous formula as in(6),

Q^{nþ1}_{S} :¼ Q^{nþ1}_{S} M W pq\ S

MðSÞ b_{q}w_{q} ð7Þ

for S = D, F and H, where Wpqrepresents the region swept out by the p-wave and q-wave over a time step Dt, and bqwqis the jump across the wave as usual.

We note that, using the same idea as described above, there is no diﬃculty to split any p-wave into m sub- waves Wpqfor q = 1, 2, . . . , m, moving oblique to the cell edge rather than normal to it. Each of these waves can then be propagated in a similar manner, computing its intersection with the neighboring cells and updat- ing the cell averages, seeFig. 1e for another illustration of a subwave, lq> 0, from the cell edge between cells C and E.

It is important to mention that the method remains fully conservative when the transverse propagation of waves is included in the algorithm. Furthermore, because waves are allowed to propagate through more than one grid cell and are not conﬁned to stay within neighboring cells, this method is typically stable as long as the time step Dt satisﬁes the usual CFL (Courant–Friedrichs–Lewy) condition relative to the regular cell sizes: Dx and Dy,

m¼ Dtmax

p;q ðkp;l_{q}Þ

minðDx; DyÞ 61; ð8Þ

even when there are very small irregular cells near the tracked interfaces. Nevertheless, the method is still only ﬁrst-order accurate, see[26,28]for more details.

2.2. High-resolution method

To achieve high resolution of the method (i.e., second-order accurate on smooth solutions, and sharp and monotone proﬁles on discontinuous solutions), as usual, we use the solution of the Riemann problem in the direction normal to each cell edge as a basis for the introduction of correction waves in a piecewise-linear form with zero mean value. Having that, we then propagate each of these waves by the corresponding speed over a time step Dt to update the cell averages it overlaps. To demonstrate this,Fig. 2a shows a typical structure of a correction wave at time tnfor a case with kp> 0 from the cell edge between regular cells G and H, i.e., the graph of the three-dimensional piecewise-linear function

Zpðx; yÞ ¼ rpðx xGÞ if ðx; yÞ 2 ½xG Dx=2; xGþ Dx=2 ½y^{a}_{GH}; y^{b}_{GH};

0 if otherwise,

ð9Þ
where rp= aprp/Dx is the slope of Zp, x_{G}is the x-coordinate of the geometric center of cell G, and y^{a}_{GH}; y^{b}_{GH}are
the lower and upper ends of the y-coordinate of this cell edge, respectively. Note that, when kp< 0, Zpcan be
made quite easily also by simply changing x_{G}in(9)to x_{H}. The new location of Zpmoving in the x-direction
over a time step Dt is Zpðx kpDt; yÞ, seeFig. 2b. The cell averages G and H are then modiﬁed by

Q^{nþ1}_{G} :¼ Q^{nþ1}_{G} 1
2

jkpjDt

Dx Dx jkpjDt
rp;
Q^{nþ1}_{H} :¼ Q^{nþ1}_{H} þ1

2 jkpjDt

Dx Dx jkpjDt rp;

ð10Þ

where the correction terms are the volume of the portion of Zpðx kpDt; yÞ (weighted by the cell size) that overlaps the grid cell. In practice, the strength of each wave should be limited by using a slope-limiter (cf.

[10,29]), and so each ap in the above rp is replaced by a limited value ea^{p} obtained by comparing ap with
the corresponding ap from the neighboring Riemann problem to the left (if kp> 0) or to the right (if
kp< 0). Conservation is clearly maintained in this step with any choice of slopes, because the above two cor-
rection terms sum to zero.

To extend the above approach to the more delicate case where the cell edge is adjacent to one or two irreg-
ular cells, we may ﬁrst deﬁne a general form of Zpin which the normal direction to the cell edge is at an arbi-
trary angle h to the x-axis. Without loss of generality, we consider the one between cells C and E shown in
Fig. 1a, as an example. Then, for the p-wave, one simple choice of the slope is rp= aprp/Dd, where Dd is some
measure of the normal distance between cells C and E, such as Dd ¼ n_{E} n_{C} with n_{S} ¼ xScos hþ y_{S}sin h rep-
resenting the n-coordinate of the geometric center of cell S in a local (n, g)-coordinate system (i.e.,
n¼ x cos h þ y sin h, g ¼ x sin h þ y cos h), for S = C, E (cf. [26]). With these deﬁnitions, for a case with
kp> 0, we form the correction wave as

Zpðn; gÞ ¼ rpn n_{0}

if ðn; gÞ 2 ½n_{0} Dd=2; n_{0}þ Dd=2 ½g^{a}_{CE};g^{b}_{CE};

0 if otherwise,

(

ð11Þ
where n0is the n-coordinate of the center of the rectangular region:½nCE;nCE Dd ½g^{a}_{CE};g^{b}_{CE}, seeFig. 2c for
an illustration. Having that, the new location of Zppropagating in the normal direction to the cell edge over a
time step Dt is Zpðn kpDt;gÞ, seeFig. 2d. Now, the cell averages aﬀected by this correction wave can be up-
dated by

Q^{nþ1}_{S} :¼ Q^{nþ1}_{S} þM Z p\ S

MðSÞ ; ð12Þ

where MðZp\ SÞ is the measure (volume) of the portion of Zpðn kpDt;gÞ that partly coincides the cell S. It is true that the unlimited slope employed in(11)would be a sensible one to use only if there is no large jump in the solution across the wave. As a general practice, it should be good also to limit the strength of each wave ap

using a slope-limiter in an upwind manner based on the sign of kp. But if we now look at the case shown in Fig. 2c, it is not obvious right away how to obtain a limited value ~ap, because we do not know exactly what is the most appropriate neighboring apto be used in the slope limiter. There are some discussions in[26]on how such corrections might be applied on irregular cells, and the algorithm would perhaps be improved by imple- menting this. However, since there is expected to be a large jump at the tracked interface, this is precisely

Fig. 2. High-resolution wave propagation method. (a) A correction wave at tnfor a case with kp> 0 from the cell edge between cells G and H. (b) Propagation of the correction wave in (a) with the speed kp> 0 over a time step Dt updates the solution in cells G and H. (c) A correction wave at tnfor a case with kp> 0 from the tracked interface between cells C and E. (d) Propagation of the correction wave in (c) with the speed kp> 0 over a time step Dt updates the solution in neighboring cells. SeeFig. 1a for the basic labeling of the grid cell.

where the limiters are expected to minimized the eﬀect of these second-order corrections and to reduce the method to Godunov. It seems reasonable to drop these corrections altogether with considerable simpliﬁcation of the algorithm.

To end this section, it should be mentioned that the wave-propagation method described here is not limited to solve hyperbolic systems in conservation form, but works equally well to solve hyperbolic systems which are not in conservation form, see[29]for the details.

3. Interface reconstruction scheme

Denote Yijas being the approximate value of the cell average of the solution(3)over the rectangular region Cij= {(x, y) | xi6x 6 xi+ Dx, yj6y 6 yj+ Dy} of grid cell (i, j). Then, when Yij2 (0, 1), we expect there is an interface (represents as a linear segment) lying somewhere within Cij. Here we are concerned with a popular piecewise linear interface construction (PLIC) method (cf.[42,45]) for the reconstruction of interfaces from the discrete set of volume fractions Y on uniform Cartesian cells which will be one of the fundamental elements in our volume tracking algorithm, see Section4. Note that in a standard PLIC method each of the linear seg- ments would be made by ﬁrst computing an approximate normal vector to the interface within the grid cell, and then determining its location in such a way that the cell volume fraction is preserved in this step. In the following, we review the aforementioned PLIC steps in more details.

3.1. Compute interface normal

As a preliminary, we consider a gradient method of Parker and Youngs[40,70]in which, when Y_{ij}2 (0, 1),
the gradient of the volume fraction at the grid cell (i, j), denoted by $Y_{ij}, is computed explicitly by a central-
diﬀerence like formula of the form

rYij¼ oY ox;oY

oy

ij

YR YL

2Dx ;YT YB

2Dy

ð13Þ with the quantities YR, YL, YT, and YB deﬁned as

Y_{R} ¼^{1}_{4}ðYiþ1;j1þ 2Yiþ1;jþ Yiþ1;jþ1Þ; Y_{L}¼^{1}_{4}ðYi1;j1þ 2Yi1;jþ Yi1;jþ1Þ;

YT ¼^{1}_{4}ðYi1;jþ1þ 2Yi;jþ1þ Yiþ1;jþ1Þ; YB¼^{1}_{4}ðYi1;j1þ 2Yi;j1þ Yiþ1;j1Þ.

Having that, the unit normal of the interface at cell (i, j) may be set to ~N_{ij}¼ rYij=jrYijj, and so in this cell

~N_{ij} points in the direction of the most rapidly decrease of Y. While the use of(13)is a very simple one for the
determination of the interface normal, it is reported in[42], for instance, that this method exhibits only overall
ﬁrst-order accuracy on the location of the interface. To achieve higher-order accuracy, such as the second
order, one popular approach is to use the least-squares method of Puckett[42,43]which is described next.

Let Eij be the sum squared error of volume fractions for the cell (i, j) of the form

Eij¼ X^{iþ1}

m¼i1

X^{jþ1}

n¼j1

Y_{mn} ~Y_{mn}

^{2}

; ð14Þ

where Ymnis the original volume fraction for cell (m, n), and ~Ymnis the approximate volume fraction for cell
(m, n) based on an extended reconstructed interface at cell (i, j). The aim of the PuckettÕs method is to obtain a
local minimum of Eij by adjusting $Yij(and so the location of the interface) iteratively in the center of a 3· 3
mesh block until convergence. Note that, during each iteration, i.e., for any given $Yij, we have to perform an
interface-ﬁnding procedure, see Section3.2, so that the condition ~Y_{ij}¼ Yij is fulﬁlled in the center cell. Here
the numerical method used to solve the minimization problem is a combination of the golden section search
and the successive parabolic interpolation (cf.[3]and the NETLIBroutineFMIN). To ensure convergence of the
correct result by the method, it is important to choose a suitable interval to work with in which Eij attains a
minimum. In practice, this is typically done by a trial-and-error procedure (cf. [15,42,50]) where an initial
guess of $Yijis picked up ﬁrst, say based on(13)or more favorably the YoungsÕ second method (see below),

and then it is expanded slowly about this initial guess to an interval until Eij at the endpoints of the interval has a bigger value than the one given by the initial guess.

Note that in PuckettÕs method $Yijis determined implicitly and also rather costly as a result of the least- squares problem(14). It is interesting to mention that a much simpler method of this type is YoungsÕ second method[71]where $Yijis obtained by solving a weighted least-squares problem of the form

E^{x}_{ij} ¼ X^{iþ1}

m¼i1

X^{jþ1}

n¼j1

xmnðYmn ~YmnÞ^{2}; ð15Þ

where in this work xmnis taken as the inverse distance from the (m, n)th cell center to the (i, j)th cell center, and
Y~_{mn}is deﬁned to be a Taylor series expansion of the volume fraction from the (m, n)th cell center to the (i, j)th
cell center, i.e.,

Y~_{mn}¼ Yijþ rYij ðxmþ1=2 xiþ1=2; y_{nþ1=2} y_{jþ1=2}Þ.

With that, the solution of the minimization problem to(15), denoted by z = ($Yij)^{T}, can be shown to satisfy
the linear system of the form

A^{T}Az¼ A^{T}b; ð16Þ

where

A¼

Dx=Ds Dy=Ds

1 0

Dx=Ds Dy=Ds

0 1

0 1

Dx=Ds Dy=Ds

1 0

Dx=Ds Dy=Ds 0

BB BB BB BB BB BB BB

@

1 CC CC CC CC CC CC CC A

; b¼

ðYi1;j1 YijÞ=Ds
ðYi1;j YijÞ=Dx
ðY_{i1;jþ1} YijÞ=Ds

ðYi;j1 YijÞ=Dy
ðYi;jþ1 YijÞ=Dy
ðY_{iþ1;j1} YijÞ=Ds

ðYiþ1;j YijÞ=Dx ðYiþ1;jþ1 YijÞ=Ds 0

BB BB BB BB BB BB BB

@

1 CC CC CC CC CC CC CC A

;

and Ds¼

ﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
ðDxÞ^{2}þ ðDyÞ^{2}
q

. In general, one may want to use either a normal equation-based or a QR-based
algorithm for the solution of this least-squares problem(15), see[63]. But after a simple algebraic calculation,
we ﬁnd that the matrix A^{T}A is simply a diagonal matrix of the form

A^{T}A¼ 2 diag 2 Dx
Ds

2

þ 1; 2 Dy Ds

2

þ 1

!

;

and so we may ﬁnd the solution of(16)very easily, and so set ~N_{ij}¼ rYij=jrYijj as the normal of the inter-
face for cell (i, j).

3.2. Determine interface location

Having gotten ~N_{ij}, we now describe a numerical procedure (cf.[45]) for an eﬃcient reconstruction of the
location of the interface over the rectangular region Cij. For convenience of the later description, let us assume
thatðx; yÞ 2 Cijis a reference point known a priori on the interface and ~Nij¼ ðaij;b_{ij}Þ with a^{2}_{ij}þ b^{2}_{ij}¼ 1. Then,
the linear-segment representation of the desired location of the interface in cell (i, j) in parametric form may be
expressed by L_{ij}ðsÞ : ðxðsÞ; yðsÞÞ ¼ ðx bijs; yþ aijsÞ, where s is a non-negative parameter value. Clearly, with
the proper choice of the interval on s, L_{ij}ðsÞ would form a line, subdividing Cijinto two parts. Denote C^{}_{ij} and
C^{þ}_{ij} to be the region occupied by subcells in the back and in the front of L_{ij}, respectively. Our aim here is to
determine ðx; yÞ and so L_{ij} in such a way that MðC^{}_{ij}Þ=MðCijÞ ¼ Yij. Surely, the latter condition ensures the
basic conservation of the volume fraction if we set Y = 1 for C^{}_{ij} and Y = 0 for C^{þ}_{ij}. Note that this setting is
consistent with our deﬁnition of the direction of ~N_{ij}, i.e., a unit vector pointing in the steepest descent direction
of Y. Recall that MðSÞ is the measure (area) of the region S.

To accomplish this, we deﬁne e^{ðkÞ}_{ij} ¼ MðC^{;k}_{ij} Þ=MðCijÞ Yijas being the error of the cell volume fraction (i, j)
at the iteration step k, where MðC^{;k}_{ij} Þ is the measure of C^{;k}_{ij} using a reference point at (x^{(k)}, y^{(k)}). With that, by
choosing a suitable initial guess: (x^{(0)}, y^{(0)}) and (x^{(1)}, y^{(1)}) such that e^{ð0Þ}_{ij} e^{ð1Þ}_{ij} <0, it is not diﬃcult to apply the
standard secant method for ﬁnding the zero of e^{ðkÞ}_{ij} , yielding the convergence of the solution
ðx^{ðkÞ}; y^{ðkÞ}Þ ! ðx; yÞ in a superlinear rate. Note that, given L_{ij}, the basic geometric information of the splitted
grid cell, such as the area and center of the cell, can be found quite easily based on algorithms in computa-
tional geometry[39].

Finally,Fig. 3shows results for a typical situation in a volume-of-ﬂuid algorithm for interface reconstruc- tion from volume fractions, where we have used both the methods of Parker and Youngs, and Puckett on a very coarse 5· 5 grid for illustrations, see [42] for an extensive study of the accuracy of these schemes and many others. Note that, in the graphs, we have set Y = 1 for grid cells in the shaded region, and Y = 0 for grid cells in the non-shadded region. As mentioned earlier, the tracked interface is approximated by a piece- wise linear segment with knots at the points where this segment intersects with the grid lines. We use the same data structure as described in[31]for the realization of the tracked interfaces in our computer program.

4. Volume tracking algorithm

In each time step, our volume-based interface tracking method for compressible ﬂuid dynamics problems in two dimensions consists of the following steps:

(1) Take a time step on the current non-uniform grid using a ﬁrst-order version of the wave-propagation method to update the cell averages of volume fractions governed by(3)at the next time step.

(2) Find the new interface location based on the volume fractions obtained in step 1 using an interface reconstruction scheme. Some cells will be subdivided, and the values in each subcell must be initialized.

(3) Take the same time interval as in step 1, but use the high-resolution wave-propagation method to update the cell averages of our model system such as(1) and (2)for single component ﬂow or a ﬂuid–mixture model proposed in[54]for multicomponent ﬂow on the new non-uniform grid created in step 2.

Note that steps 1 and 2 constitute the basic volume-moving procedure for the modeling of the propagation of interfaces, and step 3 is simply the construction of approximate solutions to the physical problems of inter- ests on a carefully chosen non-uniform grid. Since many of the key ingredients of the algorithm have been discussed before, see Sections2 and 3, our goal here is to provide some more details with application to com- pressible ﬂow problems.

4.1. Volume fraction update

Without loss of generality, suppose thatFig. 3c is the grid system we are working on at time tnwith the
approximate state for the volume fraction of (3) deﬁned as Y^{n}_{S}¼ 0 and Y^{n}_{S}¼ 1 when each of the grid cells

Fig. 3. An example of the reconstruction of the interface location from volume fractions. (a) The original cell-averaged volume fraction functions for an elliptical-shape interface. (b) Result obtained by using the method of Parker and Youngs. (c) Result obtained by using the method of Puckett. Here in the graphs we set Y = 1 and Y = 0 for grid cells in the shaded and non-shadded region, respectively.

S is in the non-shaded and shaded regions of the graph, respectively. In the meantime, the state vector for our
(single or multiple component) physical model, denoted by Q^{n}_{S}, is set by(5)to represent the approximate value
of the cell average of the solution over any given cell S. The algorithm uses a less stringent time step which is
determined by the CFL condition(8)based on the uniform grid sizes and the appropriate estimate of the wave
speeds in Q^{n}_{S}.

Notice that in this setup the value of Y^{n}is only having jump with magnitude 1 across some polygonal edges
of irregular cells, and is uniform elsewhere. Thus, in step 1, it should be reasonable to simply employ the wave-
propagation method in ﬁrst-order form (preferably the one with the transverse propagation of waves, see
Section2.1), to ﬁnd the new volume fraction Y^{n+1}on this non-uniform grid at the next time step. Here the
basic solution of the one-dimensional Riemann problem in the direction normal or tangential to the cell edge
can be computed quite easily by solving the variable-coeﬃcient linear transport equation(3)with the velocity
ﬁeld came from Q^{n} for the coeﬃcient, and Y^{n} for the initial data (cf. [29]), yielding analogously the wave
structure illustrated in Fig. 1for m = 1 over a given time step Dt.

With the above comments in mind,Fig. 4a gives numerical result of a one-step calculation for the volume fraction on the uniform rectangular grid, where the computation is performed by assuming a constant velocity ﬁeld (u, v) = (1, 1) throughout the domain, and a time step Dt = 0.06. We note that in case the grid cell is an irregular one, the result present in the ﬁgure is computed by an area-weighted average of the original data on the subcells; this is necessary in the algorithm because the reconstruction of the new interface location to be done next depends on that. It is not diﬃcult to show that our result is in good agreement as in comparison with the exact solution of this simple model problem.

4.2. New interface reconstruction

When the new volume fraction Y^{n+1}becomes available on the uniform grid, the approximate location of
the tracked interface at the next time step can be obtained in principle by applying any interface reconstruction
scheme described in Section3to the data set. If we now employ the PLIC method of Puckett to the case shown
inFig. 4a, for instance, we would have the disconnected set of linear segments as displayed inFig. 4b for the
representation of the new location of the old tracked interface over the time step Dt = 0.06. Surely, when com-
paring to the exact location of the interface, one may argue that our volume-based result does not seems to be
superior to the result gotten by using such as a more typical surface-based approach for the same problem (not
shown here, see [58] for more details). But as far as demonstrating the essential feature of the method is
concerned, this simple example serves well for our purpose here (some sample results that demonstrate the
feasibility of our method to practical problems will be present in Section5). Nevertheless, we want to empha-
size again that in our algorithm we do not assume the tracked interface to be the deﬁnitive location of the true

Fig. 4. An example of the volume-moving procedure in our interface tracking algorithm with a uniform velocity ﬁeld (u,v) = (1,1)
throughout the domain. (a) The updated volume fraction on the uniform rectangular grid based on the solution of(3)with the initial
condition shown inFig. 3c over a time step Dt = 0.06. Here j1= 5.7· 10^{3}and j2= 1.3· 10^{3}. (b) New interface location obtained by
using the PLIC method of Puckett from the volume fractions appeared in (a), where the dotted-line represents the initial location of the
tracked interface. Note that the volume fraction at the next time step is set to be unity and zero for grid cells in the shaded and non-
shadded region, respectively.

interface, but is only considered as an approximate location yielding a reﬁned grid that is better able to rep- resent the solution for the interface than the uniform grid alone.

It is clear that when a new interface is created, some cells will be subdivided, and the values in each subcell must be assigned. Here the method we use is implemented by taking a data structure that keeps the old and new grid systems separate, and as discussed in[31] doing so would be more eﬃcient than the use of a data structure for an intermediate grid that contains both the old and new subdivisions. In this manner, when the subcell is in front and back of the (directed) new tracked interface within the grid cell, we set the volume fraction at the next time step to be zero and unity, respectively. This yields the new volume fraction that will be used in the next run of our algorithm.

In addition to that, in this work we also need to initialize the state variables for the compressible model of
interests at the current time so that they can be updated in the subsequent step. For that, we should distinguish
the situation on whether the new tracked interface is inserted into a regular cell or not. Firstly, suppose the
new interface splits a regular grid cell A into two subcells B and C. In this instance, one simple way is to assign
the cell averages Q^{n}_{B} and Q^{n}_{C} based on the value Q^{n}_{A} as

Q^{n}_{B}:¼ Q^{n}_{A}; Q^{n}_{C} :¼ Q^{n}_{A}; ð17Þ

see[5,30]for some discussions on the other possibilities. Secondly, we have an irregular grid cell that is sub- divided by the old interface into subcells A and B, and also by the new interface into the subcells C and D.

Now the cell averages Q^{n}_{C} and Q^{n}_{D} may be initialized by the weighted averages of the values Q^{n}_{A} and Q^{n}_{B} as
Q^{n}_{S} :¼ ½M A \ Sð ÞQ^{n}_{A}þ M B \ Sð ÞQ^{n}_{B}=MðSÞ for S ¼ C; D. ð18Þ

4.3. Physical solution update

Our ﬁnal step of the algorithm is to ﬁnd the cell averages of the solution of our (single or multicomponent) compressible model system on the new non-uniform grid that contains some of the cell edges as the approx- imate location of the new interface at the next time step. As before (cf.[31]), this is accomplished in the same manner as described in Section2 for the solution update of hyperbolic conservation laws on a slightly non- uniform grid. That is to say, with a properly chosen Riemann solver (cf.[52–54]), we ﬁrst solve the associated one-dimensional Riemann problems in the directions normal and tangential to each of the cell edges which belongs to the original grid system when entered into the algorithm for computation. Having done so, due to the fact that we have chosen the data structure and the initializations(17) and (18)so that there are no waves generated across any polygonal edges of the irregular cell to the new grid, we may then propagate each resulting wave independently of all others to update the cell averages of the aforementioned new grid in cells neighboring to the cell edge over a given time step.

It should be remarked that as in the work of Chern and Colella[5]for interface tracking, or Berger et al.[2]

and Colella et al.[7] for handling irregular geometries in a Cartesian grid method, it is quite common to employ a post-processing step of the solution that reassigns the cell average of each irregular cells to a new one using an interpolation scheme based on the cell averages over an extended region neighboring to the grid cell under concerned. Here, in our algorithm, we take the same one-sided interpolation scheme as discussed in [5]for that matter, observing good behavior of the numerical solutions with that for many practical problems, see Section5for an example.

To end this section, we want to mention that the basic idea of the two-dimensional volume tracking algo- rithm described here can be extended in a straightforward manner to three space dimensions. That is, in steps 1 and 3 of the algorithm, we may now apply a variant of the three-dimensional wave-propagation method proposed by Langseth and LeVeque[23]on a slightly non-uniform Cartesian grid to update the volume frac- tions and the physical solutions, respectively. Clearly, one major endeavor to realize this is to devise a robust numerical procedure so that the measure (volume) of the cubical region swept out by the wave obtained from solving the one-dimensional Riemann problem over Dt and the neighboring hexahedral cells can be computed accurately and eﬃciently. It is fortunate however that in the ﬁeld of computational geometry (cf. [8,11,39]) there have been quite a few tools developed for that purpose already, and so we may simply take and imple- ment one of the popular approaches from there in our computer program. As to the interface-reconstruction

step of the algorithm, following the work appeared in[20,21], for instance, it can be done in a similar fashion without too much diﬃculty. It is important to note that to overcome the outstanding small cell problem when there is a tracked interface cutting through the underlying uniform grid, we are interested in the large time step method as usual. A detailed description of this three-dimensional algorithm (cf. [57]) will be reported else- where in the future.

5. Numerical results

We now present some sample numerical results obtained using our algorithm described in Section4for track- ing material interfaces arising from compressible multicomponent problems in two dimensions. Without stated otherwise, we have carried out all the tests in this section using the Courant number m = 0.9 deﬁned by(8), the YoungsÕ second method in the interface reconstruction scheme (for the sake of its simplicity and reasonably good prediction of the interface normal), and theMINMODlimiter in the high-resolution version of the ﬁnite- volume wave-propagation method. For comparison purposes, in many cases considered below, we have also included non-tracking results when the interface capturing version of the method is employed to the problems.

Example 5.1. To begin with, we consider a well-studied shock-bubble interaction problem that involves the
collision of a shock wave in air with a circular gas bubble (cf. [14,33,44]). To setup the problem, we take a
shock tube of size: (x, y)2 [0, 445] · [0, 89] mm^{2}, and use an initial condition that consists of a planar leftward-
moving Mach 1.22 shock wave in air located at x = 275 mm traveling from right to left, and a stationary gas
bubble with center (x0, y0) = (225, 44.5) mm and of radius r0= 25 mm in front of the shock which contains an
R22 gas inside. Here both the air and R22 are modeled as perfect gases so that we have C = c 1 (c > 1
denotes the usual ratio of speciﬁc heats) and p_{ref}= e_{ref}= 0 in the Mie–Gru¨neisen equation of state(2). For this
problem, inside the R22 gas bubble the state variables are

ðq; u; v; p; c; Y Þ ¼ ð3:863 kg=m^{3};0; 0; 1:01325 10^{5}Pa; 1:249; 1Þ;

and outside the bubble they are

ðq; u; v; p; c; Y Þ ¼ ð1:225 kg=m^{3};0; 0; 1:01325 10^{5}Pa; 1:4; 0Þ
and

ðq; u; v; p; c; Y Þ ¼ ð1:686 kg=m^{3};113:5 m=s; 0; 1:59 10^{5}Pa; 1:4; 0Þ

in the preshock and postshock regions, respectively. We note that this gives us one example that the (air–R22) material interface is accelerated by a shock wave coming from the light-ﬂuid to the heavy-ﬂuid region, and the early stage of the resulting wave pattern after the interaction would consist of a transmitted shock wave, an interface, and a reﬂected shock wave.

To solve this two-ﬂuid problem numerically, our algorithm (tracking or non-tracking version) uses a mix- ture type of the model equations proposed in [52] as the basis for approximation. Results of a sample run obtained using each of the tracking and non-tracking algorithms with an underlying uniform 3560· 356 Cartesian grid are shown inFig. 5where a sequence of the schlieren-type images of the density in some short distances around the gas bubble is presented at nine diﬀerent times t = 55, 115, 135, 187, 247, 318, 342, 417, and 1020 ls (measured relative to the time where the incident shock ﬁrst hits the upstream bubble wall). From the ﬁgure, it is easy to observe the sensible improvement of the non-tracking (i.e., interface capturing) result near the interface when the volume tracking algorithm is employed in the computation. It is also interesting to see that, as far as the qualitative structure of the solution is concerned, our tracking results are in good agree- ment with the existing experimental and numerical solutions appeared in[44, Fig. 7]. The approximate loca- tion of the tracked interfaces for the volume tracking run is displayed in Fig. 6, noticing the complicated topological changes of the interfaces during the computations. We note that it would be a notorious task to do if we want to track the interface of this degree of complexity by using an alternative surface-based inter- face tracking method like the one discussed in [31], for example.

As in[44], to get a quantitative assessment of several prominent ﬂow features, we make a diagnosis of the space–time locations of the incident shock wave, the upstream bubble wall, the downstream bubble wall, the

Fig. 5. Numerical results for a planar Mach 1.22 shock wave in air interacting with a circular R22 gas bubble. A sequence of the schlieren- type images of the density obtained using each of the tracking and non-tracking version of the methods is shown at nine diﬀerent times t = 55, 115, 135, 187, 247, 318, 342, 417, and 1020 ls with an underlying uniform 3560· 356 Cartesian grid.

refracted shock wave, and the transmitted shock wave at some selected times, and report our ﬁndings inFig. 7 only up to time t = 250 ls, seeFig. 5for an illustration of the sample positions of the marked symbols plotted in the ﬁgure and also Fig. 12 of[44]for comparison. Here we obtain the position of the incident shock, marked by symbol ‘‘·’’, by checking the jumps in the pressure to see whether there are greater than a prescribed tol- erance, at the horizontal sections of 5 mm above the bottom boundary. For the positions of the refracted shock (marked by symbol ‘‘n’’ when time t 6 202 ls), and the two transmitted shocks (marked by symbols

‘‘

’’ and ‘‘n’’ when time t > 202 ls), we do the same thing as in the incident shock case, but the pressure jumps at the axis of symmetry are checked instead. Now, as to the positions of the upstream and downstream bubble walls which are marked by symbols ’’+’’ and ‘‘e’’, respectively, there are two cases to be considered. That is, in the volume tracking run, we record the right-most and left-most of location of the tracking interface at the axis of symmetry, while in the non-tracking run we check the volume fractions to see if there are close to the midpoint value 1/2.Fig. 5 (continued)

With these trajectories, it is a common practice to perform a linear least-squares ﬁt of each set of points in Fig. 7separately, and take the slope of the respective line as one measure of the wave speed of interests. A comparison of the various computed velocities obtained using our tracking and non-tracking algorithms with those reported in the literature[14,44]is summarized inTable 1, where Vsrepresents the speed of the incident shock wave during the time t2 [0, 250] ls, Vui(Vuf) is the speed of the initial (ﬁnal) upstream bubble wall when time t2 [0, 400] ls (t 2 [400, 1000] ls), Vdi(Vdf) is the speed of the initial (ﬁnal) downstream bubble wall when time t2 [200, 400] ls (t 2 [400, 1000] ls), VRis the speed of the refracted shock when time t2 [0, 202] ls, and

time=55µs

air R22

time=115µs time=135µs

time=187µs time=247µs time=200µs

time=342µs time=417µs time=1020µs

Fig. 6. Approximate locations of the interfaces for the volume tracking run shown inFig. 5, where the dashed line in each subplot is the initial location of the air–R22 interface.

0 20 40 60

0 50 100 150 200 250

x (mm)

t (µs)

Tracking

0 20 40 60

0 50 100 150 200 250

x (mm)

t (µs)

Capturing

Fig. 7. Space-time locations of the incident shock wave (marked by symbol ‘‘·’’), the upstream bubble wall (marked by symbol ‘‘+’’), the downstream bubble wall (marked by symbol ‘‘e’’), the refracted shock (marked by symbol ‘‘n’’), and the transmitted shock (marked by symbols ‘‘’’ and ‘‘n’’) for the runs shown inFig. 5. These trajectories can be used to estimate the speed such as Vs, Vu, Vd, VR, VTin the order mentioned above.

VT is the speed of the transmitted shock when time t2 [202, 250] ls. From the table, with reference to the experimental velocities, it is clear that our algorithm (either tracking or non-tracking) gives a better prediction of the speeds for shock waves such as Vs, VR, and VT than the numerical results of Quirk and Karni. In the other instance to the prediction of the speeds for bubble walls, we still observe a better agreement of the down- stream case: Vdiand Vdf, but a slightly less accurate ones for the upstream case: Vuiand Vuf. It should be men- tioned that because the needle-shape like bubble wall in the downstream is not being tracked by the algorithm, seeFigs. 5 and 6, we do not ﬁnd a good match of the speed Vdffor that. Nevertheless, due to the sensitivity of the interface structure on the grid resolution used in the computation as well as the uncertainty on the exper- imental error, it should be fair to say that we have come up with a reasonable set of data for the bubble veloc- ities in this application.

Finally, to see how our algorithms works for the pressure, inFig. 8we record the time history of the runs shown inFig. 5at six stations, xp= 3, 11, 27, 43, 67, and 99 mm, downstream of the R22 gas bubble along the axis of symmetry. As compared this result with the one shown in Fig. 14 of[44], for example, we again notice

Table 1

A comparison of the computed velocities obtained using our tracking and non-tracking algorithms forExample 5.1with those reported in the literature[14,44]; see the text for the deﬁnition of the notations used in the table

Velocity (m/s) Vs VR VT Vui Vuf Vdi Vdf

Experiment[14] 415 240 540 73 90 78 78

Quirk and Karni[44] 420 254 560 74 90 116 82

Our result (tracking) 411 243 538 64 87 82 60

Our result (capturing) 411 244 534 65 86 98 76

100 200 300

1 2 3 4 5

xp = 3mm

p (bar)

200 250 300 350 1

1.5 2 2.5 3

xp = 11mm

200 300 400

1 1.5 2 2.5

xp = 27mm

250 300 350 400 1

1.2 1.4 1.6 1.8 2

xp = 43mm

time (µs)

p (bar)

300 400 500

1 1.2 1.4 1.6 1.8 2

xp = 67mm

time (µs) 350 400 450 500 550

1 1.2 1.4 1.6 1.8 2

xp = 99mm

time (µs)

Fig. 8. Plot of the time history of the pressure for the runs shown inFig. 5at six stations, xp= 3, 11, 27, 43, 67, and 99 mm, downstream of the R22 gas bubble along the axis of symmetry. The dotted points are the results obtained using the volume tracking method, and the solid lines are the results with interface capturing.

the sensible agreement of the solution behavior to the problem. Note that in carrying out the test we have used the solid wall boundary condition on the top and bottom, and the non-reﬂecting boundary condition on the left and right sides.

Example 5.2. Our next example of tracking material interfaces is concerned with a model underwater
explosion problem (cf.[12,55]). In this test, we take a rectangular domain (x, y)2 [2, 2] · [1.5, 1] m^{2}, and
consider the initial condition that is composed of a horizontal air–water interface at the y = 0 axis and a
circular gas bubble in water with the center (x0, y0) = (0,0.3) m and of the radius r0= 0.12 m. Here above the
air–water interface, the ﬂuid is a perfect gas at the standard atmospheric condition,

ðq; p; c; Y Þ ¼ ð1:225 kg=m^{3};1:01325 10^{5}Pa; 1:4; 1Þ.

Below the air–water interface, in region inside the gas bubble the ﬂuid is a perfect gas also with the state variables

ðq; p; c; Y Þ ¼ ð1250 kg=m^{3};10^{9}Pa; 1:4; 1Þ;

and in region outside the gas bubble the ﬂuid is water modeled by the simpliﬁed form of the Mie–Gru¨neisen equation of state(2)as

p¼ ðc 1Þqe cB; ð19Þ

where B is a pressure-like constant, and the state variables for that are
ðq; p; c; B; Y Þ ¼ ð10^{3}kg=m^{3};1:01325 10^{5}Pa; 4:4; 6 10^{8}Pa; 0Þ.

Note that, initially, all the ﬂuid components of interests is in a stationary and equilibrium state, but due to the pressure diﬀerence between the ﬂuids, breaking of the underwater bubble occurs instantaneously, yielding a circularly outward-going shock wave in water, an inward-going rarefaction wave in gas, and an interface lying in between that separates the gas and the water. Soon after this shock wave is diﬀracted through the nearby air–water interface, it is known (cf.[12,55]) that the topology of the underwater bubble will undergo a change from the original circular-shape to an oval-like shape. As time evolves, this gas bubble would con- tinue rising upward, causing the subsequent deformation of the horizontal air-water interface.

Fig. 9shows the schlieren-type images of the density and pressure at four diﬀerent times t = 0.2, 0.4, 0.8, and 1.2 ms obtained using our algorithms (tracking and non-tracking) with a 400· 250 grid. From the density plot, we clearly observe the improvement on the use of the volume tracking method over the non-tracking method to the sharpness of the solution structure near the interfaces. Moreover, from the pressure plot, we see the smooth variation of the solution near the interface, without introducing any spurious oscillations.

As far as the global structure of the solution is concerned, our result agrees quite well with the result shown in [12], for instance, where a more sophisticated surface tracking method is used in the computation. In Fig. 10, we present the location of the tracked interfaces for the volume tracking run shown in9. Note that here only the circular bubble interface is tracked initially, and the horizontal air–water interface is tracked in the algorithm only in time when its location deviates from the initial equilibrium position y = 0 axis which is aligned with the underlying grid. The cross-section of the density and pressure for the same run along line x = 0 is drawn in Fig. 11, giving some quantitative information about the diﬀerences of the tracking and non-tracking results at the selected times. The boundary conditions employed in the current computations are the solid wall boundary on the bottom, and the non-reﬂecting boundary on the remaining sides.

Example 5.3. To show how our algorithm works on tracking material interfaces with complex equation of
state, we consider the simulation of the interaction of a shock wave in molybdenum with a block of
encapsulated MORB (Mid-Ocean Ridge Basalt) liquid (cf.[37,54]). Here the computational region we take is a
unit square domain, and we use an initial condition that is composed of a planar rightward-moving Mach
1.163 shock wave in molybdenum at x = 0.3 m traveling from left to right, and a rectangular region
[0.4, 0.7]· [0, 0.5] m^{2}of MORB liquid in front of the shock. For this problem, inside the region of the MORB
liquid, we have the state variables

ðq; u; v; p; Y Þ ¼ ð2260 kg=m^{3};0; 0; 0; 1Þ;

Fig. 9. Numerical results for the simulation of a model underwater explosion problem. Schlieren-type images are shown for (a) the density, and (b) the pressure at four diﬀerent times t = 0.2, 0.4, 0.8, and 1.2 ms obtained using both the tracking and capturing version of the methods with a 400· 250 grid.

time=0.2ms

air

water

time=0.4ms

air

water

time=0.8ms

air

water

time=1.2ms

air

water

Fig. 10. Approximate locations of the tracked interface for the run shown inFig. 9, where the dashed line in each subplot is the initial location of the horizontal and circular bubble interfaces at time t = 0.

and outside the MORB, we have the state variables in the preshock region
ðq; u; v; p; Y Þ ¼ ð9961 kg=m^{3};0; 0; 0; 0Þ

and the state variables in the postshock region

ðq; u; v; p; Y Þ ¼ ð11; 042 kg=m^{3};543 m=s; 0; 3 10^{10}Pa; 0Þ.

As in[34,37,54], to model the thermodynamic behavior of solid media such as MORB and molybdenum, we use(2)together with the explicit expressions for C, pref, and erefas follows:

CðqÞ ¼ C0

q_{0}
q

a

;

p_{ref}ðqÞ ¼ p_{0}þ

q_{0}c^{2}_{0}1^{q}_{q}^{0}
1 s 1 ^{q}_{q}^{0}

h i2;

e_{ref}ðqÞ ¼ e0þ 1

2q_{0} 1q_{0}
q

p_{0}þ p_{ref}ðqÞ

ð Þ.

ð20Þ

We note that in above C0= c0 1 represents a reference Mie–Gru¨neisen coeﬃcient at q = q0, a2 [0, 1] is a dimensionless parameter, c0 denotes the zero-pressure isentropic speed of sound, and s is a dimensionless parameter which is related to the pressure derivative of the isentropic bulk modulus KS= q(op/oq)jS by (oKS/op)jS= 4s 1 (cf.[49]). With the zero initial reference state for p0and e0, typical sets of material quan- tities for the molybdenum are

ðq_{0}; c0; s;C0;aÞ ¼ ð9961 kg=m^{3};4770 m=s; 1:43; 2:56; 1Þ;

and for the MORB liquid are

ðq0; c_{0}; s;C0;aÞ ¼ ð2660 kg=m^{3};2100 m=s; 1:68; 1:18; 1Þ.

Fig. 11. Cross-sectional plots of the results for the run shown inFig. 9along x = 0, where the solid lines are the results obtained using the interface capturing method, and the dotted points are the tracking results.