Contents lists available atScienceDirect

## Journal of Computational Physics

www.elsevier.com/locate/jcp

## An Eulerian interface sharpening algorithm for compressible two-phase ﬂow: The algebraic THINC approach

### Keh-Ming Shyue

^{a}

### , Feng Xiao

^{b}

a*Institute of Applied Mathematical Sciences, National Taiwan University, Taipei 10617, Taiwan*

b*Department of Energy Sciences, Tokyo Institute of Technology, 4259 Nagatsuta Midori-ku, Yokohama 226-8502, Japan*

a r t i c l e i n f o a b s t r a c t

*Article history:*

Received 7 April 2013

Received in revised form 4 December 2013 Accepted 7 March 2014

Available online 18 March 2014

*Keywords:*

Compressible two-phase ﬂow Five-equation model Interface sharpening THINC reconstruction Mie–Grüneisen equation of state Semi-discrete wave propagation method

We describe a novel interface-sharpening approach for eﬃcient numerical resolution of a compressible homogeneous two-phase ﬂow governed by a quasi-conservative ﬁve-equation model of Allaire et al. (2001)[1]. The algorithm uses a semi-discrete wave propagation method to ﬁnd approximate solution of this model numerically. In the algorithm, in regions near the interfaces where two different ﬂuid components are present within a cell, the THINC (Tangent of Hyperbola for INterface Capturing) scheme is used as a basis for the reconstruction of a sub-grid discontinuity of volume fractions at each cell edge, and it is complemented by a homogeneous-equilibrium-consistent technique that is derived to ensure a consistent modeling of the other interpolated physical variables in the model. In regions away from the interfaces where the ﬂow is single phase, standard reconstruction scheme such as MUSCL or WENO can be used for obtaining high-order interpolated states.

These reconstructions are then used as the initial data for Riemann problems, and the resulting ﬂuctuations form the basis for the spatial discretization. Time integration of the algorithm is done by employing a strong stability-preserving Runge–Kutta method.

Numerical results are shown for sample problems with the Mie–Grüneisen equation of state for characterizing the materials of interests in both one and two space dimensions that demonstrate the feasibility of the proposed method for interface-sharpening of compressible two-phase ﬂow. To demonstrate the competitiveness of our approach, we have also included results obtained using the anti-diffusion interface sharpening method.

©2014 Elsevier Inc. All rights reserved.

**1. Introduction**

Our goal is to describe a novel Eulerian interface-sharpening approach for the eﬃcient numerical resolution of problems with material interfaces arising from inviscid compressible two-phase ﬂow. We consider an unsteady, inviscid, homogeneous two-phase ﬂow that is governed by a ﬁve-equation model system of the form

*∂*

*∂*

*t*

*(* *α*

1*ρ*

1*)* *+ ∇ · (* *α*

1*ρ*

1*u*

### *)* =

0*,*

*∂*

*∂*

*t*

*(* *α*

2*ρ*

2*)* *+ ∇ · (* *α*

2*ρ*

2*u*

### *)* =

^{0}

*,*

*E-mail addresses:*shyue@ntu.edu.tw(K.-M. Shyue),xiao@es.titech.ac.jp(F. Xiao).

http://dx.doi.org/10.1016/j.jcp.2014.03.010 0021-9991/©2014 Elsevier Inc. All rights reserved.

*∂*

*∂*

*t*

*(* *ρ*

*u*

*)* *+ ∇ · (* *ρ*

*u*

### ⊗

*u*

*)* + ∇

*p*

### =

0*,*

*∂*

*E*

*∂*

*t*

*+ ∇ · (*

^{E}*u*

### +

^{p}u### *)* =

^{0}

*,*

*∂* *α*

1
*∂*

*t*

### +

^{u}### · ∇ *α*

1### =

^{0}

*,*

(1)
as an example, for the principal motion of the state variables such as the partial densities, momentum, total energy, and volume fraction, respectively (cf.[1]). Here

*ρ*

*k*and

*α*

*k*∈ [

^{0}

*,*1]

*denote in turn the kth phasic density and volume fraction for*

*k*=

^{1}

*,*2;

*α*

1+*α*

2=^{1. We have}

*ρ*

_{=}

*α*

1*ρ*

1+*α*

2*ρ*

2 representing the total density,*u the vector of particle velocity, and p the*mixture pressure.

To close the system, the constitutive law for each of the ﬂuid phases of interest is assumed to satisfy a Mie–Grüneisen equation of state,

*p**k*

*(* *ρ*

*k*

*,*

*e*

*k*

*)* =

*p*

_{∞,}*k*

*(* *ρ*

*k*

*)* + *ρ*

*k*

*Γ*

*k*

*(* *ρ*

*k*

*)*

*e**k*

### −

*e*

_{∞,}*k*

*(* *ρ*

*k*

*)*

*,*

(2)
*where e** _{k}* is the phasic speciﬁc internal energy,

*Γ*

*k*

*= (*

^{1}

*/*

*ρ*

*k*

*)(∂p*

_{k}*/∂e*

_{k}*)|ρ*

*k*

*the Grüneisen coeﬃcient, and p*

_{∞,}

_{k}*, e*

_{∞,}*are the properly chosen states of the pressure and internal energy along some reference curve in order to match the experimental data of the material being examined. For simplicity, each of the expressions*

_{k}*Γ*

*k*

*, p*

_{∞,}*k*

*, and e*

_{∞,}*k*is taken as a function of the density only, see Section5

*for an example. As usual, E*=

*ρ*

*e*+

*ρ*

*·*

^{u}*2 is the total energy with the total internal energy deﬁned as a volume-fraction average of the form*

^{u}/*ρ*

*e*=2

*k*=^{1}

*α*

*k*

*ρ*

*k*

*e*

*k*.

*If the isobaric closure is assumed also in this model, where we have p*1=*p*2=*p in a region that contains more than*
one ﬂuid component, from the total internal energy with(2), it is easy to derive the expression for mixture pressure as

*p*

### =

### *ρ*

*e*

### −

2*k*=

^{1}

*α*

_{k}*ρ*

_{k}*e*

_{∞,}

_{k}*(* *ρ*

_{k}*)* +

2
*k*=

^{1}

*α*

_{k}

^{p}^{∞,}

^{k}^{(} *ρ*

^{(}

*k*

*)* *Γ*

_{k}*(* *ρ*

_{k}*)*

^{2}

*k*=^{1}

*α*

*k*

*Γ*

_{k}*(* *ρ*

_{k}*)* *.*

(3)
With that, it can be shown that this ﬁve-equation model is hyperbolic when each physically relevant value of the state variables of the ﬂow are deﬁned in the region of thermodynamic stability, see[1]for the detail.

As reported in the literature, conventional numerical approach originally developed for single phase compressible ﬂows can be used to solve this system. However, particular attention must be paid when computing the volume fraction. It is well known that even high order Eulerian transport scheme on ﬁxed grids cannot completely remove numerical dissipation which then tends to continuously smear out the initial jump in the volume fraction function. It is obviously problematic for immiscible interfaces. Even worse, it will eventually lead to the failure in identifying the moving interface which separate different materials for long term computations. So, we have to keep the interface sharp or at least recognizable through- out the simulation. To this end, numerical techniques, such as front tracking [5,23,45], anti-diffusion[51,47,48], interface compression[42,54], ALE (Arbitrary Lagrangian–Eulerian) and its variant[3,4,18,25,27,46,49], have been proposed.

In this work, our approach to interface-sharpening is a variant of the THINC (Tangent of Hyperbola for INterface Cap- turing) scheme that was proposed previously as an advection solver for the sharp resolution of a sub-grid discontinuity of volume fraction in the simulation of incompressible two-phase ﬂow (cf.[2,12,60–62]). Shown in these works, the improved THINC method can get numerical accuracy comparable to the VOF schemes that use explicit geometrical reconstruction, such as PLIC (piecewise linear interface calculation) reconstruction, to retrieve the moving interfaces in incompressible multi-phase ﬂows. Without the geometrical reconstruction, a THINC scheme is just a pure advection scheme, and is thus perhaps the simplest interface-capturing scheme in multi-phase ﬂow simulations.

In the present case with the ﬁve-equation model for compressible two-phase ﬂow, in regions near the interfaces, the original THINC scheme is used without any change for the reconstruction of volume fractions at each cell edge. Using the particular reconstruction function, THINC scheme effectively recovers the non-oscillatory sharp resolution of the volume fractions at every time step, and thus completely avoids the appearance of negative volume fractions that always accompa- nies the high-order Eulerian type schemes over ﬁxed computational grids. With that, assuming the retain of constant cell averages of

*ρ*

1,*ρ*

2,*u, and p within a cell, a new reconstruction procedure is derived so as to ensure a consistent modeling*of the remaining state variables such as

*α*

1*ρ*

1,*α*

2*ρ*

2,*ρ*

*u, and E at the cell edges, see*[51]for a similar procedure that was employed in an anti-diffusion type interface-sharpening algorithm. In regions away from the interfaces, where the ﬂow is single phase, standard MUSCL (Monotone Upstream-centered Schemes for Conservation Laws) or WENO (Weighted Essen- tially Non-Oscillatory) reconstruction may be applied for the interpolation of state variables to high order (cf.[20,41,57]);

this will be described further in Section4.

To ﬁnd approximate solutions of (1) with(2) and(3) numerically, we use a semi-discrete wave propagation method developed by Ketcheson et al.[15,17]for general hyperbolic systems. In this method, the spatial discretization is constructed by computing ﬂuctuations from the solutions of Riemann problems with the proposed interface-sharpened initial data at cell edges. We employ the strong stability-preserving Runge–Kutta scheme[6,7]for the integration in time, yielding a simple and eﬃcient implementation of the algorithm in the framework of the Clawpack software[16,22].

It should be mentioned that our THINC-based interface-sharpening approach is based on ﬁxed underlying grids which is different from time-varying grid approaches such as the Lagrangian or ALE[3,4,18,25,27,49], interface tracking[32,38,45],

and adaptive moving grid method[53]for sharpening interfaces. In addition, our method is simpler than the ones proposed recently by Shukla et al. [42] and So et al. [51] (see also Appendix B) in the sense that there is no need to introduce additional interface-compression and anti-diffusion source terms to the compressible two-phase ﬂow model, respectively, for the purpose of sharpening interfaces numerically (this may be diﬃcult to do for more complex multiphase ﬂow system and with phase changes, see [40,63], for example). Numerical results presented in Section 5 andAppendix A show the feasibility of the method for practical compressible two-phase ﬂow and moving interface problems.

The format of this paper is outlined as follows. In Section2, we review the basic idea of a semi-discrete ﬁnite volume method in wave-propagation form for hyperbolic conservation laws. In Section3, we describe the THINC scheme for the sharp reconstruction of interfaces based on a given set of discrete volume fractions. The interface-sharpening algorithm for compressible two-phase ﬂow is discussed in Section4. Results of some sample validation tests as well as application of the method to practical problems are shown in Section5. InAppendix A, additional advection benchmark tests are considered, and inAppendix B, a brief review of an anti-diffusion based interface sharpening scheme is described.

**2. Semi-discrete wave propagation method**

We begin our discussion by reviewing a semi-discrete ﬁnite volume method in wave-propagation form [15,17]that is essential in our interface-sharpening algorithm for solving the present quasi-conservative ﬁve-equation model numerically for compressible two-phase ﬂow. For the ease of the latter description, let us restate(1)in the form

*∂*

*q*

*∂*

*t*

### +

*N*

*j*=

^{1}

*∂*

*f*

_{j}*(*

*q*

*)*

*∂*

*x*

^{j}### +

*N*

*j*=

^{1}

*B*_{j}

*(*

*q*

*)* *∂*

*q*

*∂*

*x*

^{j}### =

^{0}

*,*

(4a)
where with*u**= ( ^{u}*1

*,u*2

*, . . . ,u*

*N*

*)the terms q, f*

*j*

*, and B*

*j*are deﬁned by

*q*

*= (* *α*

1*ρ*

1*,* *α*

2*ρ*

2*,* *ρ*

*u*1

*, . . . ,* *ρ*

*u*

*N*

*,*

*E*

*,* *α*

1*)*

^{T}*,*

(4b)
*f**j*

*= (* *α*

1*ρ*

1*u*

*j*

*,* *α*

2*ρ*

2*u*

*j*

*,* *ρ*

*u*1

*u*

*j*

### +

^{p}*δ*

*1 j*

*, . . . ,* *ρ*

*u*

*N*

*u*

*j*

### +

^{p}*δ*

*N j*

*,*

*Eu*

*j*

### +

^{pu}*j*

*,*

0*)*

^{T}*,*

(4c)
*B*_{j}

### =

^{diag}

*(*

0*,*

0*, . . . ,*

0*,*

*u*1

*δ*

_{1 j}*, . . . ,*

*u*

*N*

*δ*

_{N j}*),*

(4d)
*respectively, j*=^{1}*,*2*, . . . ,N. Hereδ**i j* *is the Kronecker delta, and N denotes the number of spatial dimension.*

*2.1. One-dimensional case*

*Consider the one-dimensional case N*=1 with a uniform grid with ﬁxed mesh spacing *x*^{1} and the number of cell in
*total M*^{1} *that discretize a spatial domain in the x*^{1} coordinate. We use a standard ﬁnite-volume formulation in which the
*value Q*_{i}*(t)approximates the cell average of the solution over the grid cell C*_{i}*: x*^{1}∈ [^{x}^{1}_{i}_{−}_{1}_{/}_{2}*,x*^{1}_{i}_{+}_{1}_{/}_{2}]^{at time t,}

*Q*_{i}

*(*

*t*

*)* ≈

^{1}

*x*

^{1}

*x*^{1}_{i}_{+}_{1}_{/}_{2}

*x*^{1}_{i}_{−}_{1}_{/}_{2}

*q*

*(*

*x*

*,*

*t*

*)*

*dx*

*.*

*Let Q(t)= ( ^{Q}*1

*,Q*2

*, . . . ,Q*

*1*

_{M}*)*

^{T}*(t)*be the vector of the approximate solution of (4)

*at time t, and*

*L*

^{1}

*(Q(t))*=

*(L*

^{1}

*(Q*1

*),L*

^{1}

*(Q*2

*), . . . ,L*

^{1}

*(Q*

*1*

_{M}*))*

^{T}*(t)*be the associated spatial-discretization vector. Then the semi-discrete version of the wave-propagation method is a method-of-lines discretization of (4) that can be written as a system of ODEs (Ordinary Differential Equations)

*∂*

*Q*

*(*

*t*

*)*

*∂*

*t*

### =

*L*

^{1}

*Q*

*(*

*t*

*)*

*,*

(5a)
where each component of*L*^{1}*(Q(t))*is deﬁned by
*L*^{1}

*Q*_{i}

*(*

*t*

*)*

### = −

^{1}

*x*

^{1}

*A*

^{+}

_{1}

*Q*

_{i}_{−}

_{1}

_{/}_{2}

### +

*A*

^{−}

_{1}

*Q*

_{i}_{+}

_{1}

_{/}_{2}

### +

*A*1

*Q*

_{i}(5b)

*for i*=^{1}*,*2*, . . . ,M*^{1}. Here*A*^{+}*Q**i*−^{1}*/*2 and *A*^{−}*Q**i*+^{1}*/*2, are the right- and left-moving ﬂuctuations, respectively, that are
entering into the grid cell, and*A*1*Q**i**is the total ﬂuctuation within C**i*. To determine these ﬂuctuations, we need to solve
Riemann problems.

Considering the ﬂuctuations *A*^{±}_{1}*Q*_{i}_{−}_{1}_{/}_{2} arising from the edge *(i*−^{1}*/*2*)between cells i*−*1 and i, for example. This*
amounts to solving a Cauchy problem that consists of

*∂*

*q*

*∂*

*t*

### + *∂*

*f*

_{1}

*(*

*q*

*)*

*∂*

*x*

^{1}

### +

*1*

^{B}*(*

*q*

*)* *∂*

*q*

*∂*

*x*

^{1}

### =

^{0}

^{(6a)}

as for the equations and

*q*

*x*

^{1}

*,*

*t*

_{0}

### =

*q*

^{L}

_{i}_{−}

_{1}

_{/}_{2}

*if x*

^{1}

*<*

*x*

^{1}

_{i}_{−}

_{1}

_{/}_{2}

*q*^{R}_{i}_{−}_{1}_{/}_{2} *if x*^{1}

*>*

*x*

^{1}

_{i}_{−}

_{1}

_{/}_{2}

*,*

^{(6b)}

*as for the initial condition at a time t*0*. Here q*^{L}_{i}_{−}_{1}_{/}_{2}=^{lim}_{x}_{→}_{x}^{1}

*(i−1/2)−** ^{q}*˜

*i*−1

*(x)and q*

_{i}

^{R}_{−}

_{1}

_{/}_{2}=

^{lim}

_{x}_{→}

_{x}^{1}

*(i−1/2)+** ^{q}*˜

*i*

*(x)*are the interpo- lated states obtained by taking limits of the reconstructed piecewise-continuous function

*q*˜

*i*−

^{1}

*(x)*or

*q*˜

*i*

*(x)*(each of them are determined based on the set of discrete data{

^{Q}*i*

*(t*

_{0}

*)}*, see Sections3and4) to the left and right of the cell edge at x

^{1}

_{i}_{−}

_{1}

_{/}_{2}, respectively.

If an approximate solver is used for the numerical resolution of the above Riemann problem (cf. [1,55]), this would
*result in three propagating discontinuities that are moving with speeds s*^{1}^{,}* ^{k}* and the jumps

*W*

^{1}

^{,}*across each of them for*

^{k}*k*=

^{1}

*,*2

*,*3, yielding the expression for the ﬂuctuations as

*A*^{±}_{1}

*Q*

_{i}_{−}

_{1}

_{/}_{2}

### =

3*k*=

^{1}

*s*

^{1}

^{,}

^{k}*q*^{L}_{i}_{−}_{1}_{/}_{2}

*,*

*q*

^{R}

_{i}_{−}

_{1}

_{/}_{2}

_{±}

*W*

^{1}

^{,}

^{k}*q*^{L}_{i}_{−}_{1}_{/}_{2}

*,*

*q*

^{R}

_{i}_{−}

_{1}

_{/}_{2}

*.*

(7a)
In a similar manner, we may deﬁne the ﬂuctuation*A*1*Q**i* *based on the Riemann problem with the initial data q*_{i}^{R}_{−}_{1}_{/}_{2} and
*q*^{L}_{i}_{+}_{1}_{/}_{2} at the cell center which gives

*A*1

*Q*

_{i}### =

3*k*=1

*s*

^{1}

^{,}

^{k}*q*^{R}_{i}_{−}_{1}_{/}_{2}

*,*

*q*

^{L}

_{i}_{+}

_{1}

_{/}_{2}

_{±}

*W*

^{1}

^{,}

^{k}*q*_{i}^{R}_{−}_{1}_{/}_{2}

*,*

*q*

_{i}

^{L}_{+}

_{1}

_{/}_{2}

### ;

^{(7b)}

this completes the deﬁnition of the ﬂuctuations and also the spatial discretization operator*L*^{1}*(Q*_{i}*(t))*in(5b). As usual, the
*notations for the quantities s*^{±}*are set by s*^{+}=^{max}*(s,*0*)and s*^{−}=^{min}*(s,*0*)*.

To integrate the system of ODEs (5a)in time, we employ the SSP (Strong Stability-Preserving) multistage Runge–Kutta scheme[6,7]. That is, in the ﬁrst-order case we use the Euler forward time discretization as

*Q*^{n}^{+}^{1}

### =

^{Q}

^{n}*+ *

^{t}L^{1}

*Q*

^{n}*,*

(8a)
*where we start with the cell average Q** ^{n}*≈

^{Q}(t

_{n}*)*

*at time t*

_{n}*, yielding the solution at the next time step Q*

^{n}^{+}

^{1}over

*t*=

*t*

*n*+

^{1}−

^{t}*n*. In the second-order case, however, we use the classical two-stage Heun method (or called the modiﬁed Euler method) as

*Q*^{∗}

### =

^{Q}

^{n}*+ *

^{t}L^{1}

*Q*

^{n}*,*

*Q*

^{n}^{+}

^{1}

### =

^{1}

2*Q*^{n}

### +

^{1}2

*Q*

^{∗}

### +

^{1}

2

*tL*

^{1}

*Q*

^{∗}

*.*

(8b)
It is common that the three-stage third-order scheme of the form

*Q*^{∗}

### =

^{Q}

^{n}*+ *

^{t}L^{1}

*Q*

^{n}*,*

*Q*

^{∗∗}

### =

^{3}

4*Q*^{n}

### +

^{1}4

*Q*

^{∗}

### +

^{1}

4

*tL*

^{1}

*Q*

^{∗}

*,*

*Q*^{n}^{+}^{1}

### =

^{1}3

*Q*

^{n}### +

^{2}

3*Q*^{∗∗}

### +

^{2}3

*tL*

^{1}

*Q*^{∗∗}

(8c) is a preferred one to be used in conjunction with the ﬁfth-order WENO scheme that is employed for the reconstruction of

˜

*q*_{i}*(x*^{1}*)*during the spatial discretization (cf.[41]).

*2.2. Multidimensional case*

To extend the one-dimensional method described above to more space dimensions, here we take a simple dimension-
by-dimension approach. Namely, assuming a uniform Cartesian grid with ﬁxed mesh spacing*x*^{1} and*x*^{2}*, in the x*^{1}-, and
*x*^{2}-direction, respectively, for example, the semi-discrete wave propagation method in two dimensions can be written of the
form

*∂*

*Q*

*(*

*t*

*)*

*∂*

*t*

### =

*L*

^{2}

*Q*

*(*

*t*

*)*

*.*

(9a)
*Here both Q(t)*= {^{Q}*i j**(t)*} ^{and}*L*^{2}*(Q(t))*= {L^{2}*(Q**i j**(t))*}*are two-dimensional arrays for i*=^{1}*,*2*, . . . ,M*^{1}*, j*=^{1}*,*2*, . . . ,M*^{2},
*with Q**i j**(t)* denoting the approximate cell average of the solution for the *(i,j)th grid cell over the region C**i j*: *(x*^{1}*,x*^{2}*)*∈
[^{x}^{1}_{i}_{−}_{1}_{/}_{2}*,x*^{1}_{i}_{+}_{1}_{/}_{2}] × [^{x}^{2}_{j}_{−}_{1}_{/}_{2}*,x*^{2}_{j}_{+}_{1}_{/}_{2}] *at time t as*

*Q*_{i j}

*(*

*t*

*)* ≈

^{1}

*x*

^{1}

*x*

^{2}

*x*^{1}_{i}_{+}_{1}_{/}_{2}

*x*^{1}_{i}_{−}_{1}_{/}_{2}
*x*^{2}_{j}_{+}_{1}_{/}_{2}

*x*^{2}_{j}_{−}_{1}_{/}_{2}

*q*

*(*

*x*

*,*

*y*

*,*

*t*

*)*

*dy dx*

*,*

and*L*^{2}*(Q**i j**(t))*representing the approximate spatial derivatives of Eqs.(4)*for C**i j* as
*L*^{2}

*Q**i j*

*(*

*t*

*)*

### = −

^{1}

*x*

^{1}

*A*

^{+}

_{1}

*Q*

*i*−

^{1}

*/*2

*,*

*j*

### +

*A*

^{−}

_{1}

*Q*

*i*+

^{1}

*/*2

*,*

*j*

### +

*A*1

*Q*

*i j*

### −

^{1}

*x*

^{2}

*A*

^{+}

_{2}

*Q*

*i*

*,*

*j*−

^{1}

*/*2

### +

*A*

^{−}

_{2}

*Q*

*i*

*,*

*j*+

^{1}

*/*2

### +

*A*2

*Q*

*i j*

### *.*

(9b)
Note that the ﬂuctuations *A*^{+}_{1}*Q**i*−^{1}*/*2*,**j*,*A*^{−}_{1}*Q**i*+^{1}*/*2*,**j* and*A*1*Q**i j* are obtained by solving the one-dimensional Rie-
*mann problems in the direction normal to the x*^{1}-axis (see (6)), while the ﬂuctuations *A*^{+}_{2}*Q*_{i}_{,}_{j}_{−}_{1}_{/}_{2}, *A*^{−}_{2}*Q*_{i}_{,}_{j}_{+}_{1}_{/}_{2} and
*A*2*Q** _{i j}* can be computed in a similar manner by solving the one-dimensional Riemann problems of in the direction nor-

*mal to the x*

^{2}-axis with

*∂*

*q*

*∂*

*t*

### + *∂*

*f*

_{2}

*(*

*q*

*)*

*∂*

*x*

^{2}

### +

*2*

^{B}*(*

*q*

*)* *∂*

*q*

*∂*

*x*

^{2}

### =

^{0}

as for the equations and with chosen reconstructed data as for the initial condition. The SSP Runge–Kutta scheme can be employed also for the integration of the semi-discrete scheme(9)in time.

It should be mentioned that the class of dimension-by-dimension semi-discrete method considered here would only give second order results when it is applied to solve general nonlinear hyperbolic systems, even if the ﬁfth-order WENO reconstruction scheme is used in the spatial discretization, and the third order Runge–Kutta method is employed in time integration, see[64]for the details. We will not discuss the fully multidimensional spatial discretization in wave-propagation form here which is beyond the scope of this paper (see[20,41]and references cited therein for more discussions on how this can be done), but will devote our attention to the devise of a new reconstruction procedure for the purpose of sharpening interfaces from compressible two-phase ﬂows.

*2.3. Include source terms*

*To end this section, we comment that if x*^{1} *is the axisymmetric direction and u*1is the radial velocity, an axisymmetric
version of the ﬁve-equation model(4)in one space dimension can be written as

*∂*

*q*

*∂*

*t*

### + *∂*

*f*1

*(*

*q*

*)*

*∂*

*x*

^{1}

### +

*B*1

*(*

*q*

*)* *∂*

*q*

*∂*

*x*

^{1}

*= ψ(*

*q*

*),*

(10)
where*ψ* is the source term derived directly from the geometric simpliﬁcation, yielding

*ψ* = − *κ*

*x*^{1}

### *α*

1*ρ*

1*u*1

*,* *α*

2*ρ*

2*u*1

*,* *ρ*

*u*

^{2}

_{1}

*,*

*Eu*1

### +

*1*

^{pu}*,*

0*T*

*.*

Here in the case of a 2D radially or 3D spherically symmetric ﬂow, the quantity

*κ*

takes 1 or 2, respectively.
Now, to ﬁnd approximate solutions of(10)numerically, in this work, we take a simple unsplit approach (cf.[20]) in that
the semi-discrete scheme(5)is employed again, but the term*L*^{1}*(Q(t))*on the right-hand side of the equation is modiﬁed
to include the effect of source term*ψ* as

*L*^{1}

*Q*

*i*

*(*

*t*

*)*

### :=

*L*

^{1}

*Q*

*i*

*(*

*t*

*)*

*+ ψ*

*Q*

*i*

*(*

*t*

*)*

*.*

(11)
We then continue by employing an SSP Runge–Kutta method as usual to integrate the resulting ODEs in time for the updated solution at the next time step, see Section5.1for a sample result obtained using the method.

**3. THINC reconstruction scheme**

The reconstruction for any physical ﬁeld to identify a sub-grid discontinuity can be derived using the THINC formulation
for volume fraction function as a building block. To describe the basic idea, we start with the one-dimensional case, and
*for simplicity rather than using x*^{1} for the spatial variable and

*α*

1*(x*

^{1}

*,t)*for the volume fraction as in Section 2, we take

*the often-employed symbols x andφ (x,t)*for that matter instead. Then as before we deﬁne the approximate value of the cell-average of

*φ (x,t)over the grid cell C*

_{i}*at time t by*

*¯φ*

*i*

*(*

*t*

*)* ≈

^{1}

*x*

*i*

*x*

*i*+1

*/*2

*x*

*i*−1

*/*2

*φ (*

*x*

*,*

*t*

*)*

*dx*

*,*

where*x**i*=^{x}*i*+^{1}*/*2−^{x}*i*−^{1}*/*2is the mesh size of the cell; assuming a non-uniform discretization of the spatial domain. Recall
that being a volume fraction of a speciﬁed ﬂuid, say ﬂuid*F*, we have the values of *¯φ**i**(t)*as follows,

*¯φ*

*i*

*(*

*t*

*)* =

### ⎧ ⎨

### ⎩

1 *if C** _{i}*is ﬁlled with ﬂuid

*F*

*,* *α* _{∈ (}

0_{∈ (}

*,*

1*)*

*if an interface lies in C*

_{i}*,*

0 if there is no ﬂuid*Fin C*

_{i}*.*

(12)

In practice, we identify an “interface cell” where the volume fraction satisﬁes,

(i)

### *< ¯* *φ*

*i*

*(*

*t*

*) <*

1### − *,*

(ii)

### ¯ *φ*

*i*+1

*(*

*t*

*)* *− ¯φ*

*i*

*(*

*t*

*) ¯* *φ*

*i*

*(*

*t*

*)* *− ¯φ*

*i*−1

*(*

*t*

*)*

*>*

0*,*

(13)
where

### is a small positive parameter (e.g., 10

^{−}

^{4}). Note that the latter condition (ii) is a monotonicity constraint on the data near the interfaces; this is introduced as a measure to prevent from incurring spurious oscillations in the method.

*Now given the volume fraction for the interface cell C**i*, we seek a reconstruction that mimics the sub-grid structure
*of the jump between 0 and 1 in the volume fraction function. For a ﬁxed time t and x*∈ [^{x}*i*−^{1}*/*2*,x*_{i}_{+}_{1}_{/}_{2}], the following
hyperbolic tangent function is well suited for this purpose,

*Φ*

*i*

*(*

*x*

*)* =

^{1}2

1

### + *σ*

*i*tanh

### *β*

*x*

### −

*x*

*i*−

^{1}

*/*2

*x*

*i*

### − ˜

^{x}*i*

*,*

(14a)
where

*σ*

*i*=

^{sgn}

*( ¯φ*

*i*

*(t))*is the sign function of the variable

*¯φ*

*i*

*(t)= ¯φ*

*i*+

^{1}

*(t)− ¯φ*

*i*−

^{1}

*(t)*, and

*β*is a prescribed parameter to control the slope and thickness of the jump which can be determined ﬂexibly in practice (here we take

*β*=

^{2}

*.*3 in all the tests considered in Section5).

Note that the only unknown in(14a)is the center* ^{x}*˜

*i*. If we assume the conservation of volume fraction of the form

1

*x*

_{i}*x*

_{i}_{+}

_{1}

_{/}_{2}

*x*_{i}_{−}_{1}_{/}_{2}

*Φ*

*i*

*(*

*x*

*)*

*dx*

*= ¯φ*

*i*

*(*

*t*

*),*

it can be determined uniquely as

### ˜

*x*

_{i}### =

^{1}

2

*β*

^{ln}exp

*(β(*

1### + *σ*

*i*

### −

^{2}

*¯φ*

_{i}*)/* *σ*

*i*

*)*

1### −

^{exp}

*(β(*

1### − *σ*

*i*

### −

^{2}

*¯φ*

_{i}*)/* *σ*

*i*

*)*

*.*

(14b)
It should be mentioned that using (14) as a building block, we may reconstruct any physical ﬁeld *ψ(x,t)* *for x*∈
[^{x}*i*−^{1}*/*2*,x**i*+^{1}*/*2]where it has a jump*ψ**i**= ψ*max*− ψ*minby,

*Ψ*

*i*

*(*

*x*

*)* *= ψ*

min*+ Φ*

*i*

*(*

*x*

*)ψ*

*i*

*.*

(15)
Here lower and upper bounds of the jump*ψ*min and*ψ*max can be determined from the neighboring cells such as*ψ*min=
min*(ψ*_{i}^{L}_{−}_{1}_{/}_{2}*, ψ*_{i}^{R}_{+}_{1}_{/}_{2}*)* and*ψ*max=^{max}*(ψ*_{i}^{L}_{−}_{1}_{/}_{2}*, ψ*_{i}^{R}_{+}_{1}_{/}_{2}*)*, where *ψ*_{i}^{L}_{−}_{1}_{/}_{2} denotes the value computed from the reconstruction
over the left-side cell[^{x}*i*−^{3}*/*2*,x**i*−^{1}*/*2]^{and}*ψ*_{i}^{R}_{+}_{1}_{/}_{2} from the right-side cell[^{x}*i*+^{1}*/*2*,x**i*+^{3}*/*2]^{.}

Reconstruction(15)can be applied to various physical variables, such as primitive variables, characteristic variables and
ﬂux functions. As long as the THINC reconstruction(15)is built up, we can get directly the values at the cell boundaries as
*a function of time t through either a point mapping or a trajectory integration for a wave propagation (transport) equation.*

In the latter instance, we consider the following wave equation as an example,

*∂ψ*

*∂*

*t*

### +

^{u}*∂ψ*

*∂*

*x*

### =

^{0}

*,*

*where u is the characteristic velocity ﬁeld. We split u into right-moving and left-moving parts, i.e., u*=^{u}^{+}+^{u}^{−} ^{where}
*u*^{+}=^{max}*(u,*0*)and u*^{−}=^{min}*(u,*0*). The solution at the right-side of cell boundary x*_{i}_{−}_{1}_{/}_{2}reads

*ψ*

_{i}

^{R}_{−}

_{1}

_{/}_{2}

*(*

*t*

*)* *= Ψ*

*i*

*x*

_{i}_{−}

_{1}

_{/}_{2}

### −

^{u}^{−}

_{i}_{−}

_{1}

_{/}_{2}

^{t}*= ψ*

min*+ Φ*

*i*

*x*

*i*−

^{1}

*/*2

### −

*u*

^{−}

_{i}_{−}

_{1}

_{/}_{2}

*t*

*ψ*

*i*

*= ψ*

min### +

^{1}2

1

### + *σ*

*i*tanh

### *β*

*−*

_{u}*i*−

^{1}

*/*2

*t*

*x*

*i*

### − ˜

^{x}*i*

*ψ*

*i*

*,*

(16a)
*while that at the left-side of cell boundary x**i*+^{1}*/*2reads

*ψ*

_{i}

^{L}_{+}

_{1}

_{/}_{2}

*(*

*t*

*)* *= Ψ*

*i*

*x*

_{i}_{+}

_{1}

_{/}_{2}

### −

^{u}^{+}

_{i}_{+}

_{1}

_{/}_{2}

^{t}*= ψ*

min*+ Φ*

*i*

*x*

_{i}_{+}

_{1}

_{/}_{2}

### −

^{u}^{+}

_{i}_{+}

_{1}

_{/}_{2}

^{t}*ψ*

_{i}*= ψ*

min### +

^{1}2

1

### + *σ*

*i*tanh

### *β*

1### −

^{u}+*i*+^{1}*/*2*t*

*x*

_{i}### − ˜

*x*

*i*

*ψ*

*i*

*,*

(16b)
which can be used to any approximate Riemann solver in a full- or semi-discrete form, see the formulas given in[2]for that purpose also.

For the mass conservation law,

*∂ψ*

*∂*

*t*

### + *∂*

*f*

*∂*

*x*

### =

0*,*

*where f* =^{u}ψ*is the ﬂux function, we have the total leftward ﬂux across cell boundary x*_{i}_{−}_{1}_{/}_{2}*over t as,*

### ˆ

*f*

_{i}^{−}

_{−}

_{1}

_{/}_{2}

*(*

*t*

*)* =

*t*0

*f*

*(*

*x*

_{i}_{−}

_{1}

_{/}_{2}

*,*

*t*

*)*

*dt*

### =

*x*

*i*−

^{1}

*/*2

*x*

_{i}_{−}

_{1}

_{/}_{2}−

^{u}^{−}

_{i}_{−}

_{1}

_{/}_{2}

^{t}*Ψ*

_{i}*(*

*x*

*)*

*dx*

*= ψ*

0*u*

^{−}

_{i}_{−}

_{1}

_{/}_{2}

*t*

### − *σ*

*i*

*x*

*2*

_{i}*β*

^{ln}

_{cosh}

*(β(˜*

^{x}*i*

### +

^{u}^{−}

^{i}

_{}^{−}

^{1}

_{x}

^{/}

_{i}^{2}

^{t}*))*

cosh*(β* ˜

*x*

_{i}*)*

*ψ*

_{i}*.*

(17a)
*Analogously, the rightward ﬂux across x**i*+^{1}*/*2is

### ˆ

*f*

_{i}^{+}

_{+}

_{1}

_{/}_{2}

*(*

*t*

*)* =

*t*0

*f*

*(*

*x*

*i*+

^{1}

*/*2

*,*

*t*

*)*

*dt*

### =

*x*

*i*+1

*/*2

*x*

_{i}_{+}

_{1}

_{/}_{2}−

^{u}^{+}

_{i}_{+}

_{1}

_{/}_{2}

^{t}*Ψ*

*i*

*(*

*x*

*)*

*dx*

*= ψ*

0*u*

^{+}

_{i}_{+}

_{1}

_{/}_{2}

*t*

### − *σ*

_{i}*x*

*2*

_{i}*β*

^{ln}

_{cosh}

*(β(*

1### − ˜

^{x}*i*

### −

^{u}^{+}

^{i}

_{}^{+}

^{1}

_{x}

^{/}

_{i}^{2}

^{t}*))*

cosh*(β(*

1### − ˜

^{x}*i*

*))*

*ψ*

_{i}*.*

(17b)
where*ψ*0*= (ψ*min*+ ψ*max*)/*2. So, a one-step ﬁnite volume formulation to transport mass can be written as,

*¯ψ*

_{i}*(*

*t*

*)* *= ¯ψ*

*i*

*(*

0*)* −

^{1}

*x*

_{i}### ˆ

_{f}

_{i}_{+}

_{1}

_{/}_{2}

*(*

*t*

*)* − ˆ

^{f}*i*−1

*/*2

*(*

*t*

*)* *,*

which is used in the existing THNIC method for incompressible ﬂows.

To extend the above THINC reconstruction scheme to more than one space dimension, here we take a simple dimension- by-dimension approach as proposed in Xiao et al.[60]; this one-dimensional version of the method works quite well with the semi-discrete wave propagation method described in Section2and also the interface-sharpening reconstruction scheme for compressible two-phase ﬂow described next.

**4. Interface sharpening algorithm**

In each time step, our interface-sharpening algorithm for compressible two-phase ﬂow problems consists of the following steps:

(1) Reconstruct a piecewise polynomial function, denoted by* ^{q}*˜

*(x*

^{i}*,t*

_{n}*), for all x*

*based on standard MUSCL/WENO recon-*

^{i}*struction procedure from the cell average Q*

^{n}*at time t*

*to more than ﬁrst order.*

_{n}(2) Modify*q*˜*(x*^{i}*,t**n**)for interface cells using a variant of THINC scheme from Q** ^{n}* to sharpen the resolution of interfaces.

(3) Solve Riemann problems with interpolated initial data from*q*˜*(x*^{i}*,t**n**)*obtained in steps 1 and 2 for spatial discretization.

*(4) Employ a semi-discrete method in wave propagation form to update Q*^{n}*from the current time to the next Q*^{n}^{+}^{1} over
a time step*t.*

Note that if step 2 is omitted in the above algorithm, it is simply the standard semi-discrete method proposed by Ketche- son et al.[15,17]for general hyperbolic systems with new applications to compressible two-phase ﬂow, see Section2. In this instance, the implementation of the MUSCL/WENO reconstruction scheme as suggested in[13,20,34]for step 1 is enough for that matter. Our goal here is to describe step 2 of the reconstruction procedure that is essential to our interface-sharpening algorithm.

*4.1. Homogeneous-equilibrium-consistent reconstruction scheme*

We begin step 2 by employing the one-dimensional THINC scheme described in Section3for the reconstruction of the
*volume fraction function q ^{(}*

^{N}^{+}

^{4}

*=*

^{)}*α*

1, in the ﬁve-equation model(4). Suppose that with Q

^{n}*given at a time t*

*n*, we have

*obtained the interpolated volume fraction for the ith cell at the edges(*

*α*

1*)*

_{i}

^{R}_{−}

_{1}

_{/}_{2}and

*(*

*α*

1*)*

_{i}

^{L}_{+}

_{1}

_{/}_{2}based on the relations(16a) and(16b)

*at t*=0, respectively.

With that, to construct the sub-grid structure of the partial density, we follow the approach proposed by So et al.[51]

(see alsoAppendix Band [54]) in that the phasic densities

*ρ*

1 and*ρ*

2 are assumed to remain constants within the cells
during the reconstruction step, yielding the deﬁnition of the cell-edge states for the partial densities, i.e., the ﬁrst and
*second component of q; q*

^{(}

^{k}*=*

^{)}*α*

_{k}*ρ*

_{k}*for k*=

^{1}

*,*2, as

*(* *α*

*k*

*ρ*

*k*

*)*

_{i}

^{R}_{−}

_{1}

_{/}_{2}

*= (* *α*

*k*

*ρ*

*k*

*)*

*i*

*+ (* *ρ*

*k*

*)*

*i*

### *(* *α*

*k*

*)*

_{i}

^{R}_{−}

_{1}

_{/}_{2}

*− (* *α*

*k*

*)*

*i*

### *,*

(18a)
*(* *α*

_{k}*ρ*

_{k}*)*

^{L}

_{i}_{+}

_{1}

_{/}_{2}

*= (* *α*

_{k}*ρ*

_{k}*)*

*i*

*+ (* *ρ*

_{k}*)*

*i*

### *(* *α*

_{k}*)*

_{i}

^{L}_{+}

_{1}

_{/}_{2}

*− (* *α*

_{k}*)*

*i*

### *.*

(18b)
Here*(*

*α*

*k*

*)*

*i*

*(*

*ρ*

*k*

*)*

*i*, and

*(*

*α*

*k*

*ρ*

*k*

*)*

*i*

*are the cell-average variables obtained from Q*

_{i}

^{n}*in an interface cell C*

*.*

_{i}*Now to ﬁnd the reconstructed states for the total momentum q ^{(}*

^{j}^{+}

^{2}

*=*

^{)}*ρ*

*u*

_{j}*for j*=

^{1}

*,*2

*, . . . ,N, the velocity ﬁeld*

*=*

^{u}*(u*1

*,u*2

*, . . . ,u*

*N*

*)*is assumed to remain unchanged within the cells also; this is true in a region for the contact discontinuity (the case we are interested in this work). Thus, we have

*(* *ρ*

*u*

*)*

^{R}

_{i}_{−}

_{1}

_{/}_{2}

*= (* *ρ*

*u*

*)*

_{i}### +

^{u}*i*

### *ρ*

_{i}

^{R}_{−}

_{1}

_{/}_{2}

_{−} *ρ*

_{i}^{} *,*

(18c)
*(* *ρ*

^{u}*)*

^{L}

_{i}_{+}

_{1}

_{/}_{2}

*= (* *ρ*

^{u}*)*

*i*

### +

^{u}*i*

### *ρ*

_{i}

^{L}_{+}

_{1}

_{/}_{2}

### − *ρ*

*i*

### *,*

(18d)
where

*ρ*

_{i}

^{R}_{−}

_{1}

_{/}_{2}

_{=}

^{}

^{2}

_{k}_{=}

_{1}

*(*

*α*

*k*

*ρ*

*k*

*)*

^{R}

_{i}_{−}

_{1}

_{/}_{2}and

*ρ*

_{i}

^{L}_{+}

_{1}

_{/}_{2}

_{=}

^{}

^{2}

_{k}_{=}

_{1}

*(*

*α*

*k*

*ρ*

*k*

*)*

^{L}

_{i}_{+}

_{1}

_{/}_{2}denotes the reconstructed mixture density at the left- and

*right-edge of the grid cell i, respectively.*

*Finally, to reconstruct the total energy q ^{(}*

^{N}^{+}

^{3}

*=*

^{)}*E, we assume further the pressure equilibrium within interface cells,*yielding

*E*^{R}_{i}_{−}_{1}_{/}_{2}

### =

^{E}*i*

### +

^{K}*i*

### *ρ*

_{i}

^{R}_{−}

_{1}

_{/}_{2}

### − *ρ*

*i*

### +

2*k*=1

*(* *ρ*

*k*

*e*

_{k}*)*

*i*

### *(* *α*

*k*

*)*

^{R}

_{i}_{−}

_{1}

_{/}_{2}

*− (* *α*

*k*

*)*

*i*

### *,*

(18e)
*E*^{L}_{i}_{+}_{1}_{/}_{2}

### =

^{E}*i*

### +

^{K}*i*

### *ρ*

_{i}

^{L}_{+}

_{1}

_{/}_{2}

### − *ρ*

*i*

### +

2*k*=1

*(* *ρ*

*k*

*e*

_{k}*)*

*i*

### *(* *α*

*k*

*)*

^{L}

_{i}_{+}

_{1}

_{/}_{2}

*− (* *α*

*k*

*)*

*i*

### *,*

(18f)
*where K** _{i}*=

^{u}*i*·

^{u}*i*

*/2 denotes the speciﬁc kinetic energy in cell C*

*. Note that since the above sub-grid reconstruction proce- dure(18)makes use of the basic homogeneous-equilibrium assumptions of the underlying continuum model(1), it should be pertinent to call it a homogeneous-equilibrium-consistent reconstruction scheme.*

_{i}*4.2. Interface only problem*

To see how our algorithm works for sharpening interfaces, it is instructive to consider an interface only problem where
*the initial data consists of uniform pressure p*=* ^{p}*0, constant velocity

*u*=

*0, and constant phasic densities*

^{u}*ρ*

_{k}_{=}

*ρ*

*for*

_{k0}*k*=

^{1}

*,*2 throughout the domain, while there are jumps on the other variables such as partial densities and volume fractions

*across some interfaces. Without loss of generality, we consider a one-dimensional problem with a positive velocity u*

_{1}=

*u*0

*>*0. In this case, suppose that a 3-wave HLLC approximate solver (cf. [9,44,55]) is used in step 3 of the algorithm for solving the Riemann problem (6) (similar result follows if a Roe solver [1] is employed instead), it is easy to ﬁnd the

*ﬂuctuations for each cell C*

_{i}*, i*=

^{1}

*,*2

*, . . . ,M*

^{1}, as

*A*^{−}_{1}

*Q*

*i*+

^{1}

*/*2

### =

0*,*

*A*

^{+}

_{1}

*Q*

*i*−

^{1}

*/*2

### =

*u*0

*q*

_{i}

^{R}_{−}

_{1}

_{/}_{2}

### −

*q*

^{L}

_{i}_{−}

_{1}

_{/}_{2}

### *,*

*A*1

*Q*

_{i}### =

*0*

^{u}*q*

^{L}

_{i}_{+}

_{1}

_{/}_{2}

### −

^{q}

^{R}

_{i}_{−}

_{1}

_{/}_{2}

### *.*

Inserting the above expression for ﬂuctuations into(5b), we have the spatial discretization of this interface only problem:

*L*^{1}

*Q*

_{i}*(*

*t*

*)*

### = −

^{1}

*x*

^{1}

*u*

_{0}

*q*_{i}^{L}_{+}_{1}_{/}_{2}

### −

^{q}

_{i}

^{L}_{−}

_{1}

_{/}_{2}

### *.*

Now if the Euler method(8a)is employed for the time integration in(5a), together with the above spatial discretization
*term the cell average Q*_{i}* ^{n}*is updated by

*Q*^{n}_{i}^{+}^{1}

### =

^{Q}

_{i}

^{n}### −

*t*

*x*

^{1}

*u*

_{0}

*q*_{i}^{L}_{+}_{1}_{/}_{2}

### −

^{q}

_{i}

^{L}_{−}

_{1}

_{/}_{2}