An Eulerian interface sharpening algorithm for compressible two-phase flow: The algebraic THINC approach

Contents lists available atScienceDirect

Journal of Computational Physics

www.elsevier.com/locate/jcp

An Eulerian interface sharpening algorithm for compressible two-phase flow: The algebraic THINC approach

Keh-Ming Shyue

^a

, Feng Xiao

^b

aInstitute of Applied Mathematical Sciences, National Taiwan University, Taipei 10617, Taiwan

bDepartment 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 flow 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 efficient numerical resolution of a compressible homogeneous two-phase flow governed by a quasi-conservative five-equation model of Allaire et al. (2001)[1]. The algorithm uses a semi-discrete wave propagation method to find approximate solution of this model numerically. In the algorithm, in regions near the interfaces where two different fluid 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 flow 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 fluctuations 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 flow. 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 efficient numerical resolution of problems with material interfaces arising from inviscid compressible two-phase flow. We consider an unsteady, inviscid, homogeneous two-phase flow that is governed by a five-equation model system of the form

∂

∂

t

( α

1

ρ

1

) + ∇ · ( α

1

ρ

1u

 ) =

0

,

∂

∂

t

( α

2

ρ

2

) + ∇ · ( α

2

ρ

2u

 ) =

⁰

,

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

+

^pu

 ) =

⁰

,

∂ α

1

∂

t

+ 

^u

· ∇ α

1

=

⁰

,

(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∈ [⁰,1]denote in turn the kth phasic density and volume fraction for k=¹,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 fluid phases of interest is assumed to satisfy a Mie–Grüneisen equation of state,

pk

( ρ

k

,

ek

) =

p_∞,k

( ρ

k

) + ρ

k

Γ

k

( ρ

k

) 

ek

−

e_∞,k

( ρ

k

) 

,

(2)

where e_k is the phasic specific internal energy,Γk= (¹/

ρ

k)(∂p_k/∂e_k)|ρk the Grüneisen coefficient, and p_∞,k, e_∞,k 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, p_∞,k, and e_∞,kis taken as a function of the density only, see Section5 for an example. As usual, E=

ρ

e+

ρ

^u· ^u/2 is the total energy with the total internal energy defined as a volume-fraction average of the form

ρ

e=2

k=¹

α

k

ρ

kek.

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

p

=

 ρ

e

−



2 k=¹

α

_k

ρ

_ke_∞,k

( ρ

_k

) +



2 k=¹

α

_k^p∞,k

⁽ ρ

k

) Γ

_k

( ρ

_k

)

 

²

k=¹

α

k

Γ

_k

( ρ

_k

) .

(3)

With that, it can be shown that this five-equation model is hyperbolic when each physically relevant value of the state variables of the flow are defined 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 flows 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 fixed 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 flow (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 flows. 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 flow simulations.

In the present case with the five-equation model for compressible two-phase flow, 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 fixed 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 flow 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 find 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 fluctuations 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 efficient 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 fixed 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 flow model, respectively, for the purpose of sharpening interfaces numerically (this may be difficult to do for more complex multiphase flow 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 flow and moving interface problems.

The format of this paper is outlined as follows. In Section2, we review the basic idea of a semi-discrete finite 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 flow 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 finite volume method in wave-propagation form [15,17]that is essential in our interface-sharpening algorithm for solving the present quasi-conservative five-equation model numerically for compressible two-phase flow. For the ease of the latter description, let us restate(1)in the form

∂

q

∂

t

+



N j=¹

∂

f_j

(

q

)

∂

x^j

+



N j=¹

B_j

(

q

) ∂

q

∂

x^j

=

⁰

,

(4a)

where withu= (^u1,u2, . . . ,uN)the terms q, fj, and Bjare defined by

q

= ( α

1

ρ

1

, α

2

ρ

2

, ρ

u1

, . . . , ρ

uN

,

E

, α

1

)

^T

,

(4b)

fj

= ( α

1

ρ

1uj

, α

2

ρ

2uj

, ρ

u1uj

+

^p

δ

1 j

, . . . , ρ

uNuj

+

^p

δ

N j

,

Euj

+

^puj

,

0

)

^T

,

(4c)

B_j

=

^diag

(

0

,

0

, . . . ,

0

,

u1

δ

_{1 j}

, . . . ,

uN

δ

_{N j}

),

(4d)

respectively, j=¹,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 fixed mesh spacing x¹ and the number of cell in total M¹ that discretize a spatial domain in the x¹ coordinate. We use a standard finite-volume formulation in which the value Q_i(t)approximates the cell average of the solution over the grid cell C_i: x¹∈ [^x1_i−1/2,x¹_i+1/2]^{at time t,}

Q_i

(

t

) ≈

¹



x¹

x¹_i+1/2

x¹_i−1/2

q

(

x

,

t

)

dx

.

Let Q(t)= (^Q1,Q2, . . . ,Q_M1)^T(t) be the vector of the approximate solution of (4) at time t, and L¹(Q(t))= (L¹(Q1),L¹(Q2), . . . ,L¹(Q_M1))^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¹



Q

(

t

) 

,

(5a)

where each component ofL¹(Q(t))is defined by L¹



Q_i

(

t

) 

= −

¹



x¹



A⁺₁



Q_i−1/2

+

A⁻₁



Q_i+1/2

+

A1



Q_i



(5b)

for i=¹,2, . . . ,M¹. HereA⁺Qi−¹/2 and A⁻Qi+¹/2, are the right- and left-moving fluctuations, respectively, that are entering into the grid cell, andA1Qiis the total fluctuation within Ci. To determine these fluctuations, we need to solve Riemann problems.

Considering the fluctuations A^±₁Q_i−1/2 arising from the edge (i−¹/2)between cells i−1 and i, for example. This amounts to solving a Cauchy problem that consists of

∂

q

∂

t

+ ∂

f₁

(

q

)

∂

x¹

+

^B1

(

q

) ∂

q

∂

x¹

=

^{0 (6a)}

as for the equations and

q



x¹

,

t₀



=

q^L_i−1/2 if x¹

<

x¹_i−1/2

q^R_i−1/2 if x¹

>

x¹_i−1/2

,

^(6b)

as for the initial condition at a time t0. Here q^L_i−1/2=^lim_x→x¹

(i−1/2)−^q˜i−1(x)and q_i^R_−1/2=^lim_x→x¹

(i−1/2)+^q˜i(x)are the interpo- lated states obtained by taking limits of the reconstructed piecewise-continuous functionq˜i−¹(x)orq˜i(x)(each of them are determined based on the set of discrete data{^Qi(t₀)}, see Sections3and4) to the left and right of the cell edge at x¹_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 jumpsW^1,k across each of them for k=¹,2,3, yielding the expression for the fluctuations as

A^±₁



Q_i−1/2

=



3 k=¹



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 define the fluctuationA1Qi 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

A1



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 definition of the fluctuations and also the spatial discretization operatorL¹(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 first-order case we use the Euler forward time discretization as

Qⁿ⁺¹

=

^Qn

+ 

^tL¹



Qⁿ



,

(8a)

where we start with the cell average Qⁿ≈^Q(t_n) at time t_n, yielding the solution at the next time step Qⁿ⁺¹ over t= tn+¹−^tn. In the second-order case, however, we use the classical two-stage Heun method (or called the modified Euler method) as

Q^∗

=

^Qn

+ 

^tL¹



Qⁿ



,

Qⁿ⁺¹

=

¹

2Qⁿ

+

¹ 2Q^∗

+

¹

2



tL¹



Q^∗



.

(8b)

It is common that the three-stage third-order scheme of the form

Q^∗

=

^Qn

+ 

^tL¹



Qⁿ



,

Q^∗∗

=

³

4Qⁿ

+

¹ 4Q^∗

+

¹

4



tL¹



Q^∗



,

Qⁿ⁺¹

=

¹ 3Qⁿ

+

²

3Q^∗∗

+

² 3



tL¹



Q^∗∗



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

˜

q_i(x¹)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 fixed mesh spacingx¹ andx², in the x¹-, and x²-direction, respectively, for example, the semi-discrete wave propagation method in two dimensions can be written of the form

∂

Q

(

t

)

∂

t

=

L²



Q

(

t

) 

.

(9a)

Here both Q(t)= {^Qi j(t)} ^andL²(Q(t))= {L²(Qi j(t))}are two-dimensional arrays for i=¹,2, . . . ,M¹, j=¹,2, . . . ,M², with Qi j(t) denoting the approximate cell average of the solution for the (i,j)th grid cell over the region Ci j: (x¹,x²)∈ [^x1_i−1/2,x¹_i+1/2] × [^x2_j−1/2,x²_j+1/2] at time t as

Q_{i j}

(

t

) ≈

¹



x¹



x²

x¹_i+1/2

x¹_i−1/2 x²_j+1/2

x²_j−1/2

q

(

x

,

y

,

t

)

dy dx

,

andL²(Qi j(t))representing the approximate spatial derivatives of Eqs.(4)for Ci j as L²



Qi j

(

t

) 

= −

¹



x¹



A⁺₁



Qi−¹/2,j

+

A⁻₁



Qi+¹/2,j

+

A1



Qi j



−

¹



x²



A⁺₂



Qi,j−¹/2

+

A⁻₂



Qi,j+¹/2

+

A2



Qi j

 .

(9b)

Note that the fluctuations A⁺₁Qi−¹/2,j,A⁻₁Qi+¹/2,j andA1Qi j are obtained by solving the one-dimensional Rie- mann problems in the direction normal to the x¹-axis (see (6)), while the fluctuations A⁺₂Q_i,j−1/2, A⁻₂Q_i,j+1/2 and A2Q_{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²-axis with

∂

q

∂

t

+ ∂

f₂

(

q

)

∂

x²

+

^B2

(

q

) ∂

q

∂

x²

=

⁰

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 fifth-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 flows.

2.3. Include source terms

To end this section, we comment that if x¹ is the axisymmetric direction and u1is the radial velocity, an axisymmetric version of the five-equation model(4)in one space dimension can be written as

∂

q

∂

t

+ ∂

f1

(

q

)

∂

x¹

+

B1

(

q

) ∂

q

∂

x¹

= ψ(

q

),

(10)

whereψ is the source term derived directly from the geometric simplification, yielding

ψ = − κ

x¹

 α

1

ρ

1u1

, α

2

ρ

2u1

, ρ

u²₁

,

Eu1

+

^pu1

,

0



T

.

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

κ

takes 1 or 2, respectively.

Now, to find 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 termL¹(Q(t))on the right-hand side of the equation is modified to include the effect of source termψ as

L¹



Qi

(

t

) 

:=

L¹



Qi

(

t

) 

+ ψ 

Qi

(

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 field 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¹ for the spatial variable and

α

1(x¹,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 define the approximate value of the cell-average ofφ (x,t)over the grid cell C_i at time t by

¯φ

i

(

t

) ≈

¹



xi x

i+1/2 xi−1/2

φ (

x

,

t

)

dx

,

wherexi=^xi+¹/2−^xi−¹/2is the mesh size of the cell; assuming a non-uniform discretization of the spatial domain. Recall that being a volume fraction of a specified fluid, say fluidF, we have the values of ¯φi(t)as follows,

¯φ

i

(

t

) =

⎧ ⎨

⎩

1 if C_iis filled with fluidF

, α _{∈ (}

0

,

1

)

if an interface lies in C_i

,

0 if there is no fluidFin C_i

.

(12)

In practice, we identify an “interface cell” where the volume fraction satisfies,

(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⁻⁴). 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 Ci, we seek a reconstruction that mimics the sub-grid structure of the jump between 0 and 1 in the volume fraction function. For a fixed time t and x∈ [^xi−¹/2,x_i+1/2], the following hyperbolic tangent function is well suited for this purpose,

Φ

i

(

x

) =

¹ 2



1

+ σ

itanh

 β



x

−

xi−¹/2



xi

− ˜

^xi



,

(14a)

where

σ

i=^sgn( ¯φi(t))is the sign function of the variable ¯φi(t)= ¯φi+¹(t)− ¯φi−¹(t), andβ is a prescribed parameter to control the slope and thickness of the jump which can be determined flexibly in practice (here we takeβ=².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

=

¹

2

β

^ln



exp

(β(

1

+ σ

i

−

²

¯φ

_i

)/ σ

i

)

1

−

^exp

(β(

1

− σ

i

−

²

¯φ

_i

)/ σ

i

)



.

(14b)

It should be mentioned that using (14) as a building block, we may reconstruct any physical field ψ(x,t) for x∈ [^xi−¹/2,xi+¹/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[^xi−³/2,xi−¹/2]^andψ_i^R_+1/2 from the right-side cell[^xi+¹/2,xi+³/2]^.

Reconstruction(15)can be applied to various physical variables, such as primitive variables, characteristic variables and flux 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

=

⁰

,

where u is the characteristic velocity field. 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/2reads

ψ

_i^R_−1/2

(

t

) = Ψ

i



x_i−1/2

−

^u−_i−1/2^t



= ψ

min

+ Φ

i



xi−¹/2

−

u⁻_i−1/2t



ψ

i

= ψ

min

+

¹ 2



1

+ σ

itanh

 β



_u− i−¹/2t



xi

− ˜

^xi



ψ

i

,

(16a)

while that at the left-side of cell boundary xi+¹/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

+

¹ 2



1

+ σ

itanh

 β



1

−

^u

+i+¹/2t



x_i

− ˜

xi



ψ

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 flux function, we have the total leftward flux across cell boundary x_i−1/2over t as,

ˆ

f_i⁻_−1/2

(

t

) =

t 0

f

(

x_i−1/2

,

t

)

dt

=

x

i−¹/2 x_i−1/2−^u−_i−1/2^t

Ψ

_i

(

x

)

dx

= ψ

0u⁻_i−1/2t

− σ

i



x_i 2

β

^ln



_cosh

(β(˜

^xi

+

^u−i_⁻¹_x^/_i^2t

))

cosh

(β ˜

x_i

)



ψ

_i

.

(17a)

Analogously, the rightward flux across xi+¹/2is

ˆ

f_i⁺_+1/2

(

t

) =

t 0

f

(

xi+¹/2

,

t

)

dt

=

x

i+1/2 x_i+1/2−^u+_i+1/2^t

Ψ

i

(

x

)

dx

= ψ

0u⁺_i+1/2t

− σ

_i



x_i 2

β

^ln



_cosh

(β(

1

− ˜

^xi

−

^u+i_⁺¹_x^/_i^2t

))

cosh

(β(

1

− ˜

^xi

))



ψ

_i

.

(17b)

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

¯ψ

_i

(

t

) = ¯ψ

i

(

0

) −

¹



x_i

 ˆ

_fi+1/2

(

t

) − ˆ

^fi−1/2

(

t

)  ,

which is used in the existing THNIC method for incompressible flows.

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 flow described next.

4. Interface sharpening algorithm

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

(1) Reconstruct a piecewise polynomial function, denoted by^q˜(xⁱ,t_n), for all xⁱ based on standard MUSCL/WENO recon- struction procedure from the cell average Qⁿat time t_n to more than first order.

(2) Modifyq˜(xⁱ,tn)for interface cells using a variant of THINC scheme from Qⁿ to sharpen the resolution of interfaces.

(3) Solve Riemann problems with interpolated initial data fromq˜(xⁱ,tn)obtained in steps 1 and 2 for spatial discretization.

(4) Employ a semi-discrete method in wave propagation form to update Qⁿ from the current time to the next Qⁿ⁺¹ over a time stept.

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 flow, 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 five-equation model(4). Suppose that with Qⁿ given at a time tn, 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 definition of the cell-edge states for the partial densities, i.e., the first and second component of q; q^(k)=

α

_k

ρ

_kfor k=¹,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)iare the cell-average variables obtained from Q_iⁿ in an interface cell C_i.

Now to find the reconstructed states for the total momentum q^(j+2)=

ρ

u_j for j=¹,2, . . . ,N, the velocity field ^u= (u1,u2, . . . ,uN)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

+ 

^ui

 ρ

_i^R_−1/2

₋ ρ

_i

^ ,

(18c)

( ρ 

^u

)

^L_i+1/2

= ( ρ 

^u

)

i

+ 

^ui

 ρ

_i^L_+1/2

− ρ

i

 ,

(18d)

where

ρ

_i^R_−1/2=^2_k=1(

α

k

ρ

k)^R_i−1/2and

ρ

_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

=

^Ei

+

^Ki

 ρ

_i^R_−1/2

− ρ

i

 +



2 k=1

( ρ

ke_k

)

i

 ( α

k

)

^R_i−1/2

− ( α

k

)

i

 ,

(18e)

E^L_i+1/2

=

^Ei

+

^Ki

 ρ

_i^L_+1/2

− ρ

i

 +



2 k=1

( ρ

ke_k

)

i

 ( α

k

)

^L_i+1/2

− ( α

k

)

i

 ,

(18f)

where K_i= ^ui· ^ui/2 denotes the specific kinetic energy in cell C_i. 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.

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=^p0, constant velocity u= ^u0, and constant phasic densities

ρ

_k=

ρ

_k0 for k=¹,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₁= u0>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 find the fluctuations for each cell C_i, i=¹,2, . . . ,M¹, as

A⁻₁



Qi+¹/2

=

0

,

A⁺₁



Qi−¹/2

=

u0



q_i^R_−1/2

−

q^L_i−1/2

 ,

A1



Q_i

=

^u0



q^L_i+1/2

−

^qR_i−1/2

 .

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

L¹



Q_i

(

t

) 

= −

¹



x¹u₀



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ⁿis updated by

Qⁿ_i⁺¹

=

^Q_iⁿ

− 

t



x¹u₀



q_i^L_+1/2

−

^q_i^L_−1/2



,