A modified HLLC-type Riemann solver for the compressible six-equation two-fluid model

Computers & Fluids 76 (2013) 86–104

Contents lists available at SciVerse ScienceDirect

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

A modified HLLC-type Riemann solver for the compressible six-equation two-fluid model Geum-Su Yeom a,⇑, Keun-Shik Chang b a

School of Mechanical & Automotive Engineering, Kunsan National University, Miryong-dong, Gunsan-si, Jeonbuk 573-701, South Korea Division of Aerospace Engineering, Department of Mechanical, Aerospace and Systems Engineering, Korea Advanced Institute of Science and Technology, 373-1 Gusong-dong, Yuseung-gu, Daejeon-si 305-701, South Korea

b

a r t i c l e

i n f o

Article history: Received 11 October 2010 Received in revised form 26 December 2012 Accepted 31 January 2013 Available online 16 February 2013 Keywords: Six-equation two-fluid model HLLC scheme HLL scheme Two-phase shock tube Shock–bubble interaction

a b s t r a c t We present a modified HLLC-type Riemann solver that is applicable to the compressible six-equation two-fluid model of the two-phase flow. The modified HLLC scheme restores all the characteristic fields that have been neglected by the recent HLLC scheme proposed by Zein et al. (Zein A, Hantke M, Warnecke G. Modeling phase transition for compressible two-phase flows applied to metastable liquids. J Comput Phys 2010;229:2964–98). We test the modified HLLC scheme on several 1D two-phase shock tube problems and a 2D shock–bubble interaction problem. The modified HLLC scheme is proved more accurate than other schemes like the Zein et al.’s HLLC scheme and the earlier two-phase HLL scheme developed by the present authors (Yeom GS, Chang KS. Numerical simulation of two-fluid two-phase flows by HLL scheme using an approximate Jacobian matrix. Numer Heat Transfer B-Fund 2006;49:155–77). The numerical test indicates that the present HLLC scheme is very robust and effective in acquiring high resolution for the waves in the six-equation two-fluid model. Ó 2013 Elsevier Ltd. All rights reserved.

1. Introduction Numerical simulation of two-phase flow is important for both physical scientists and engineers because the two-phase flows are frequently encountered in nature and technology. Wide applications can be found in such industry as petroleum processing, chemical and mechanical engineering, and nuclear power plant. Nevertheless, development of mathematical modeling and numerical methods on the two-phase flow has never been sufficient because of some inherent difficulties. There have appeared a variety of models for simulation of the two-phase flow: see Ishii and Hibiki [1], Saurel and Lemetayer [2] for general description. Among them, we focus on the two-fluid models that are derived by an ensemble average of the local instantaneous phasic flow equations. In general, the six-equation two-fluid model [3–9] and the seven-equation two-fluid model [2,10–12] have been well accepted by the two-phase flow research community. The six-equation model consists of two sets of the three 1D conservation laws that govern the physics of the liquid and the gas phases. Here, an assumption is made that the pressures of the two fluids are locally in equilibrium resulting in a single pressure ⇑ Corresponding author. Tel.: +82 64 469 4712, mobile: +82 10 3474 5082; fax: +82 64 469 4727. E-mail address: [email protected] (G.-S. Yeom). 0045-7930/$ - see front matter Ó 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.compfluid.2013.01.021

and two velocities. In contrast, the seven-equation model includes the void fraction evolution equation in addition to the above six conservation laws. This model introduces the concept of two pressures and two velocities. In order to meet the interface conditions on the boundary of the two fluids, relaxation operators for pressure and velocity are necessary [2,11]. Furthermore, a thermal relaxation [12] can be used when heat exchange between phases is taken into account. The seven-equation model is easier and more adapted to the upwind numerical schemes since this model enjoys the analytic eigenvalues and eigenvectors of the hyperbolic equation system [2,11]. In contrast, the six-equation model is computationally more efficient because it dispenses with the additional void fraction evolution equation. Further, the six-equation model can skip the pressure relaxation since it inherently assumes a single equilibrium pressure. In this paper we formulate an accurate HLLC-type method for the six-equation two-fluid model. The original thermo-hydraulic system codes such as RELAP5, CATHARE, and TRAC are based on the six-equation two-fluid model, using the donor-cell differencing method on the staggered grids. Unfortunately, the equation system as well as the scheme has showed poor numerical stability and resolution for the strong thermal and mechanical non-equilibrium problems [3]. The virtual mass and the interfacial pressure terms have been added to the equations for stability [3,6] that availed the more accurate approximate Riemann solvers (see Toro [13]). For example, Städtke et al.

87

G.-S. Yeom, K.-S. Chang / Computers & Fluids 76 (2013) 86–104

Fig. 1. Wave structures for the six-equation two-fluid model.

[4] implemented the Flux Vector Splitting method to solve the transient two-phase flows. Sainsaulieu [14], Toumi [3], and Cortes [15] used the Roe-type method by using approximate eigenvalues and eigenvectors. They derived the generalized Rankine–Hugoniot conditions for the non-conservative system. Coquel et al. [16] successfully proposed a two-step approach by decoupling the hydrodynamic system and keeping equality of pressures. The AUSM-family scheme was successfully implemented to the seven-equation two-fluid model by Niu [34]. He solved dilute particulate turbulence flows through a 90° bend. Subsequently, Paillère et al. [17] and Niu [5] developed the AUSM-family schemes for the six-equation two-fluid model. The HLL scheme was first presented by Harten et al. [18] to approximately solve the Riemann problem for the hyperbolic equation system. The central idea of the HLL scheme is to approximate the speeds of the two fastest characteristic waves in 1D while the intermediate wave (contact discontinuity) is neglected. Saurel et al. [2,11] applied the HLL scheme to the seven-equation two-fluid model. More recently, Yeom and Chang [6–9] formulated the HLL scheme for the six-equation two-fluid model to solve a few 1D as well as 2D problems. The scheme has proved robust and efficient despite some inherent diffusion because of the intermediate waves ignored in the model. An improved HLLC scheme was introduced by Toro et al. [19] for the Euler equations that restored the intermediate wave neglected by the HLL scheme. Recent development shows an extension of the HLLC scheme to the two-fluid models. Saurel et al. [24] implemented the HLLC scheme to the single velocity, pressure non-equilibrium model (the so-called six-equation model) that is derived from the seven-equation model of Baer and Nunziato [10] in the asymptotic limit of stiff velocity relaxation only. Tokareva and Toro [20] developed the HLLC-type Riemann solver for the Baer–Nunziato model [10] of compressible two-phase flows, considering all the characteristic waves in the exact Riemann solution. Unfortunately, the scheme requires an expensive Newton–Raphson iteration to determine the unknown variables of the intermediate states. Further, a HLLC-type scheme has been developed for the sevenequation two-fluid model by Zein et al. [12] who has accounted the intermediate wave in their numerical model and used the pressure, velocity, and thermal relaxation operators. Nonetheless, the number of waves in their 1D model was only three, resembling the HLLC scheme applied to the Euler equations, not seven for the sevenequation two-fluid equation model. Therefore, the Zein et al.’s HLLC scheme still ignores some characteristic fields in the two fluids. To the best of the authors’ knowledge, no HLLC-type scheme has been completely formulated for the six-equation two-fluid model. In this paper, we develop a modified HLLC-type Riemann solver for the six-equation two-fluid model by taking account of all the characteristic fields in the system. For this purpose, we first introduce the idea of the recent HLLC-type scheme by Zein et al. [12] to

Fig. 2. Comparison of the wave structure between the Zein et al.’s HLLC-type Riemann solver and the modified HLLC-type Riemann solver for the two-fluid model. (a) Zein’s HLLC-type Riemann solver. (b) Present HLLC-type Riemann solver.

the present six-equation two-fluid model. Then we improve the scheme by restoring the intermediate waves neglected in the earlier work. We then solve one-dimensional two-phase shock tube problems and a two-dimensional shock–bubble interaction problem to show the accuracy and the robustness of the proposed scheme. 2. Six-equation two-fluid model The compressible six-equation two-fluid two-phase flow model can be written as follows:

@U @F @ ag þ þH ¼ S; @t @x @x

ð1Þ

with

1 0 1 0 ag qg ug ag qg 1 U1 B U C B ag q ug C B a q u2 þ a p C g g g g g C C B 2C B B g C C B C B B B U 3 C B ag qg Eg C B ag qg ug Eg þ ag ug pg C C C C; B B B ; F ¼ U¼B C¼B C C B al ql ul C B U 4 C B a l ql C B C C B C B B @ U 5 A @ al ql ul A @ al ql u2l þ al pl A 0

al ql El U6 1 1 0 0 0 B pi C B FD C C C B B B i iC B i D C B p u C B uF C C C B B H¼B C; S ¼ B 0 C: C B 0 C B C B i C B @ p A @ F D A pi ui ui F D 0

al ql ul El þ al ul pl

We use notations as follows. a: void fraction, q: density, u: velocity, E: specific total energy, p: pressure, pi: interfacial pressure, ui: interfacial velocity, and FD: interfacial friction. The subscript g is used for the gas and l for the liquid. In the present work the phase change and heat transfer are not considered. The specific total energy E is defined by

88

G.-S. Yeom, K.-S. Chang / Computers & Fluids 76 (2013) 86–104

Fig. 3. Abgrall’s test problem with a stationary contact discontinuity. 2

E¼eþ

u ; 2

ð2Þ

ð4Þ

The interfacial friction term is given in the general form as

where e is the specific internal energy. The void fractions satisfy

ag þ al ¼ 1:

pg ¼ pl ¼ p:

ð3Þ

The gas and liquid pressures are assumed equal in this model,

F D ¼ C D ag al qg ðug  ul Þ;

ð5Þ

where a positive constant CD is the drag coefficient. In this paper, we assume the drag coefficient is very large, i.e., CD  1, so that the gas

G.-S. Yeom, K.-S. Chang / Computers & Fluids 76 (2013) 86–104

89

Fig. 4. Abgrall’s test problem with a moving contact discontinuity.

and liquid velocities become equal very quickly. In this case, we can apply the instantaneous velocity relaxation proposed by Saurel et al. [2,11] for integrating the interfacial friction term.

system. As a remedy, the virtual mass or the interfacial pressure terms would be included in the governing equations [3,6]. In this paper, we use instead the interfacial pressure model [21,6] defined by

2.1. Interfacial pressure and velocities

pi ¼ p  d

As stated earlier, the six-equation two-fluid model could result in numerical instability because it does not guarantee a hyperbolic

where d = 1.2 is the damping coefficient. Note that the equation system becomes hyperbolic for d P 1 if the relative velocity is small.

ag al qg ql ðu  ul Þ2 ; ag ql þ al qg g

ð6Þ

90

G.-S. Yeom, K.-S. Chang / Computers & Fluids 76 (2013) 86–104

Fig. 5. Cavitation problem: cavitation region at center is developed by opposite expansion waves.

The interfacial velocity ui is applied at the center of mass and modeled after Saurel and Abgrall [11],

ag qg ug þ al ql ul u ¼ : ag qg þ al ql i

ð7Þ

Each phase is treated as compressible fluid, governed by the modified stiffened-gas equations of state (EOS) [29].

ð8Þ ð9Þ

where k = 1, 2. Cv,k is the heat capacity at constant volume. ck, p1,k, and qk are the EOS parameters. Note that the above EOS becomes the conventional stiffened-gas EOS when qk = 0. The speed of sound is defined by

ck ¼

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ck ðpk þ p1;k Þ

qk

:

From the conservative variable vector U (=Ui, i = 1, 2, . . . , 6), we get the six primitive variables (a, p, u, q, E and e) for the two fluids. The velocities are given by

ug ¼

2.2. Equation of state

pk þ ck p1;k ek ¼ þ qk ; ðck  1Þqk pk þ p1;k ; Tk ¼ C v ;k qk ðck  1Þ

2.3. Primitive variables from the conservative variables

ð10Þ

U2 ; U1

ul ¼

U5 : U4

ð11Þ

The specific total energies are

Eg ¼

U3 ; U1

El ¼

U6 : U4

ð12Þ

The specific internal energies are

eg ¼ Eg 

u2g ; 2

e l ¼ El 

u2l : 2

ð13Þ

We derive the following relation from the void fraction relation, Eq. (3), the equal pressure assumption, Eq. (4), and the EOS, Eq. (8):

1 ¼ ag þ al ¼

¼

U1

qg

þ

U4

ql

U 1 ðcg  1Þðeg  qg Þ U 4 ðcl  1Þðel  ql Þ þ : p þ cl p1;l p þ cg p1;g

ð14Þ

G.-S. Yeom, K.-S. Chang / Computers & Fluids 76 (2013) 86–104

91

Fig. 6. The high-pressure water/low-pressure air shock tube problem with a pressure ratio of 104.

Eq. (14) gives a second-order polynomial equation with respect to the common pressure p. By taking its positive root, we get the analytic form of the pressure as

1h p ¼ eg U 1 þ qg U 1  cg ðp1;g  eg U 1 þ qg U 1 Þ  el U 4 2 pffiffiffiffii þql U 4  cl ðp1;l  el U 4 þ ql U 4 Þ þ D ;

ð15Þ

where

D ¼ ½cg p1;g þ cl p1;l þ ðcg  1Þðeg  qg ÞU 1 2 þ 2ðcl  1Þðel  ql Þ

ð16Þ Unþ1 ¼ Uni  i

Finally, the void fractions are given by

qg

ð20Þ

The homogenous system, Eq. (18), is discretized using the finite volume method as

Using the EOS, we can compute the densities as

U1

ð19Þ

3.2. Hyperbolic operator

 ql Þ2 U 24 :

ag ¼

ð18Þ

where LH is the hyperbolic operator used in the homogenous Eq. (18), and LS is the source operator for the interfacial friction source Eq. (19). The source operator can be solved by the conventional Runge–Kutta method or the instantaneous velocity relaxation procedure [2,11] when the drag coefficient is very large. The solution for n + 1 time step is then obtained by

Unþ1 ¼ LS LH Un :

 ½cg p1;g  cl p1;l þ ðcg  1Þðeg  qg ÞU 1 U 4 þ ðcl  1Þ2 ðel

qg ¼ qg ðp; eg Þ; ql ¼ ql ðp; el Þ:

@U @F @ ag þ þH ¼ 0; @t @x @x @U ¼ S; LS : @t

LH :

 Dt h i Dt   Fiþ1=2  Fi1=2  Hni ðag Þiþ1=2  ðag Þi1=2 ; Dx Dx ð21Þ

;

al ¼ 1  ag :

ð17Þ

3. Numerical method 3.1. Fractional step method We use the fractional step method, decomposing Eq. (1) to the two subsystems as

Fiþ1=2

 g Þiþ1=2

where and ða are respectively the numerical flux and the numerical void fraction at the cell interface which can be calculated by some approximate Riemann solvers. In order to use the HLLC-type Riemann solver, the eigenstructure of the equation system needs to be determined a priori. Unfortunately, its analytic eigenvalues have not been obtained for the six-equation model that was made hyperbolic by including the interfacial pressure term. Instead, the approximate but analytic eigenvalues have been evaluated in different forms by several

92

G.-S. Yeom, K.-S. Chang / Computers & Fluids 76 (2013) 86–104

Fig. 7. Effect of a small volume fraction e for the water–air shock tube problem.

authors [3,4,6,7,14,22]. For instance, Yeom and Chang [6] derived the analytic eigenvalues based on the approximate Jacobian matrix as

k1 ¼ u g ; k2 ¼ u l ;

able to assume that the two linearly degenerate characteristic fields is very close to each other, i.e., Sg⁄  Sl⁄  S⁄ because ug  ul. In this section, we first implement the Zein’s HLLC-type Riemann solver to the present six-equation two-fluid model. The numerical flux at the cell interface then becomes

k3 ¼ u g  c g ; k4 ¼ u g þ c g ; sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ag cg c2l p1;l ; k5 ¼ u l  c l 1  ag cg cl ðcl  1Þp1;l þ al qg c2g cl þ ag ql c2l cg sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ag cg c2l p1;l k6 ¼ u l þ c l 1  ; ag cg cl ðcl  1Þp1;l þ al qg c2g cl þ ag ql c2l cg where c is the ratio of specific heats and p1 is the parameter of the stiffened-gas equation of state. It is known that the two k1- and k2- characteristic fields represented by Sg⁄ and Sl⁄ in Fig. 1 are the linearly degenerate ones and the rest characteristic fields, SlL, SgL, SgR, and SlR, are genuinely nonlinear [22,6]. 3.2.1. Zein’s HLLC-type Riemann solver Recently, Zein et al. [12] proposed the HLLC-type Riemann solver for the two-fluid model with the velocity, pressure, and thermal relaxation operators by considering a wave structure consisting of three characteristic fields, i.e., SL, S⁄, and SR: see Fig. 2a. The Zein’s HLLC Riemann solver [12] is based on the following assumption: when the drag force between phases is very large so that the two phasic velocities are almost identical, it is reason-

Fiþ1=2

8 FL > > > < FL þ SL ðUL  UL Þ ¼ > F > R þ SR ðUR  UR Þ > : FR

for 0 6 SL ; for SL 6 0 6 S ; for S 6 0 6 SR ;

ð22Þ

for 0 P SR :

The intermediate states, U⁄L and U⁄R, are

0

agK qgK ðSK  ugK Þ

1

C B SK  S C B C B a q ðS  u Þ gK gK K gK C B S C B  SK S B " !# C C B C B agK qgK ðSK  ugK Þ pgK C B E þ ðS  u Þ S þ gK  gK  B SK  S qgK ðSK  ugK Þ C C B UK ¼ B C; C B a q ðS  u Þ lK lK K lK C B C B S  S K  C B C B alK qlK ðSK  ulK Þ C B C B S C B S  S K  C B     A @ alK q ðSK  ulK Þ plK lK ElK þ ðS  ulK Þ S þ SK  S qlK ðSK  ulK Þ ð23Þ where the index K = L or R.

93

G.-S. Yeom, K.-S. Chang / Computers & Fluids 76 (2013) 86–104

Fig. 8. The effect of the EOS parameters for the water–air shock tube problem.

The numerical void fraction is given by  g Þiþ1=2

ða



agL for 0 6 S ; ¼ agR for 0 P S :

ð24Þ

Then, the left and the right fastest wave speeds can be approximated by

SL ¼ minðugL  cgL ; ulL  clL ; ugR  cgR ; ulR  clR Þ; SR ¼ maxðugL þ cgL ; ulL þ clL ; ugR þ cgR ; ulR þ clR Þ:

ð25Þ ð26Þ

The speed of the intermediate wave S⁄ is calculated by

S ¼

pR  pL þ qL uL ðSL  uL Þ  qR uR ðSR  uR Þ ; qL ðSL  uL Þ  qR ðSR  uR Þ

Table 1 Summary of the stiffened-gas EOS parameters for liquid water used by earlier researchers. Case

Authors

cl

p1,l

1

Saurel and Abgrall [11] Saurel and LeMetayer [2]

4.4

6  108

2

Paillère et al. [17] Liou et al. [32]

2.8

8.5  108

3

Chang and Liou [30] Niu et al. [5]

1.932

1.645  109

ð27Þ

In the above, we use the mixture values for the pressure, velocity, and density variables as below:

pK ¼ agK pgK þ alK plK ;

ð28Þ

qK ¼ agK qgK þ alK qlK ; agK qgK ugK þ alK qlK ulK uK ¼ ; qK

ð29Þ ð30Þ

where K = L or R. 3.2.2. A modified HLLC-type Riemann solver Though the Zein’s HLLC-type Riemann solver resolves the linearly degenerate characteristic field (contact discontinuity) accurately, the model incorporates only three waves instead of the

full seven waves only: one intermediate wave for the linearly degenerate filed plus two fastest waves for the genuinely nonlinear characteristic fields. This wave model is clearly insufficient since the problem has many other nonlinear waves induced by the two-fluid model. It is comparable to the HLL scheme attempting to resolve the three wave system by two fastest waves for the 1D Euler equations. In this paper, we present an improved version of the Zein’s HLLC-type Riemann solver [12] that restores all the missing waves for the genuinely nonlinear characteristic fields. As indicated in the earlier papers [12,20,23], it is an important property that the void fraction changes only across the contact wave S⁄. On this reason, the two fluids can be decoupled across the wave S⁄ and the two-fluid governing equations are reduced to a pair of two single-fluid equations (Euler equations) [20]. Based

94

G.-S. Yeom, K.-S. Chang / Computers & Fluids 76 (2013) 86–104

Fig. 9. The high-pressure water/low-pressure air shock tube problem with a pressure ratio of 105.

on this observation, we reinstate all the missing intermediate states by considering complete wave structures for the six-equation two-fluid model. Fig. 2b shows a wave configuration of the modified HLLC-type Riemann solver, where the conservative variables of the two fluids are denoted by Ug and Ul as

0

ag qg

1

B C Ug ¼ @ ag qg ug A; ag qg Eg

0

al ql

0 B

UgK ¼

agK qgK ðSgK  ugK Þ B B SgK  S

¼

ðFg Þiþ1=2 ðFl Þiþ1=2

0

ð31Þ UlK ¼

ð32Þ

;

alK qlK ðSlK  ulK Þ B B

ðFl Þiþ1=2

8 FlL > > > < F þ S ðU  U Þ lL lL lL lL ¼ > FlR þ SlR ðUlR  UlR Þ > > : FlR

pgK

gK ðSgK ugK Þ

for 0 6 SgL ; for SgL 6 0 6 S ; for S 6 0 6 SgR ; for 0 P SgR ;

1

S h @ ElK þ ðS  ulK Þ S þ q

plK lK ðSlK ulK Þ

C C i A;

ð36Þ



agL for 0 6 S ; agR for 0 P S :

ð37Þ

The speeds of the four nonlinear waves are estimated as

ð33Þ

for 0 6 SlL ; for SlL 6 0 6 S ; for S 6 0 6 SlR ;

SlK  S

1

where K = L or R. The numerical void fraction at the cell interface is given by

ðag Þiþ1=2 ¼

where

ðFg Þiþ1=2

C C C; C iA ð35Þ

!

8 FgL > > > < F þ S ðU  U Þ gL gL gL gL ¼ > F gR þ SgR ðUgR  UgR Þ > > : FgR

EgK

S h þ ðS  ugK Þ S þ q

1

B C Ul ¼ @ al ql ul A: a l q l El

The resulting numerical flux of the modified HLLC-type Riemann solver is given by

Fiþ1=2

B @

1

1

ð34Þ

for 0 P SlR :

The intermediate states, Ug⁄L, Ug⁄R, Ul⁄L, and Ul⁄R are given by

SgL ¼ minðugL  cgL ; ugR  cgR Þ;

ð38Þ

SgR ¼ maxðugL þ cgL ; ugR þ cgR Þ;

ð39Þ

SlL ¼ minðulL  clL ; ulR  clR Þ;

ð40Þ

SlL ¼ maxðulL þ clL ; ulR þ clR Þ:

ð41Þ

The speed of the contact wave S⁄ is computed by a simple average of two contact wave speeds for each fluid as

S ¼

1 ðSg þ Sl Þ; 2

ð42Þ

G.-S. Yeom, K.-S. Chang / Computers & Fluids 76 (2013) 86–104 Table 2 Parameters of the modified stiffened-gas equation of state for vapor and liquid dodecanes. Phase

c

p1 (Pa)

Cv (J/kg K)

Cp (J/kg K)

q (J/kg)

Vapor Liquid

1.025 2.35

0 4  108

1.956  103 1.077  103

2.005  103 2.535  103

237  103 755  103

95

where the individual contact wave speed is estimated using the mixture variables in a way comparable to Zein et al. [12], namely,

pR  pL þ qL uL ðSgL  uL Þ  qR uR ðSgR  uR Þ ; qL ðSgL  uL Þ  qR ðSgR  uR Þ p  pL þ qL uL ðSlL  uL Þ  qR uR ðSlR  uR Þ ; Sl ¼ R qL ðSlL  uL Þ  qR ðSlR  uR Þ

Sg ¼

ð43Þ ð44Þ

Fig. 10. Comparison of the present results with Zein et al.’s results for the Zein’s dodecane liquid–vapor shock tube problem. The results are obtained with 1250 cells. (a) Zein et al.’s results [12]: the Baer–Nunziato type model. (b) Present results: the six-equation model.

96

G.-S. Yeom, K.-S. Chang / Computers & Fluids 76 (2013) 86–104

Fig. 11. Comparison of the present HLLC scheme with Zein’s HLLC scheme for the Zein’s dodecane liquid–vapor shock tube problem, solved by the six-equation two-fluid model.

where the mixture values are defined by Eqs. (28)–(30) again. Note that the present approximation, S ¼ 12 ðSg þ Sl Þ, was taken due to its simplicity, but the other approximations such as Eq. (27) by Zein et al. [12] are possible. According to a wide range of our numerical tests, the present approximation agrees very well with the exact solution and is very robust. The two-dimensional extension of the present scheme can be easily made by the finite volume method with the numerical flux and the numerical void fraction evaluated at the cell interface. For instance, see Appendix D of Saurel et al. [24] or Yeom and Chang [7]. The second-order extension both in time and space is made by the weighted average flux method [25,9].

ab ¼ ai ; qb ¼ qi ; pb ¼ pi ; ub ¼ ui þ 2up ; v b ¼ v i :

Here, subscripts i and b represent the interior and the ghost cells, respectively. The time step Dt is determined by

Dt ¼ CFL 

Dx ; kmax

The boundary conditions are imposed on the ghost cells that are the mirror image of the boundary cells. The numerical flux at the boundary is then computed using the ghost cell values by the same Riemann solver used for the internal cells. For the transmissive boundaries, the ghost cell values are taken the same as the boundary cell values,

ab ¼ ai ; qb ¼ qi ; pb ¼ pi ; ub ¼ ui ; v b ¼ v i :

ð45Þ

For the piston boundary with a prescribed speed up, the ghost cell values are given as

ð47Þ

where CFL 2 (0, 1) is the Courant–Friedrichs–Lewy number. The maximum wave speed at the current time level, kmax, is computed from

kmax ¼ maxðjug j þ cg ; jul j þ cl Þ: 3.3. Boundary and CFL condition

ð46Þ

ð48Þ

The actual time step is then taken with the minimum among all the cells,

Dt ¼ mini ðDti Þ for i ¼ 1; 2; . . . ; I:

ð49Þ

4. Numerical tests 4.1. Abgrall’s test problems We first consider a void fraction jump in the two uniform twophase flows: one is the stationary contact problem and the other is

97

G.-S. Yeom, K.-S. Chang / Computers & Fluids 76 (2013) 86–104

Fig. 12. Effect of e for the Zein’s dodecane liquid–vapor shock tube problem.

the moving contact problem. These problems are designed to check agreement with the Abgrall’s principle [26]: A flow, uniform in pressure and velocity, must remain uniform in the same variables during its time evolution. A horizontal pipe, 10 m long, has a diaphragm at 3 m separating the left and the right initial states given below:

ag 1 0 0:2 B p C B 105 Pa B C B B C B B ug C B 10 m=s B C ¼B Bu C B 10 m=s B lC B B C B @ qg A @ 1:3 kg=m3 0

ql 1

0

0

ag 0:2 B p C B 105 Pa B C B B C B B ug C B 0 m=s B C ¼B Bu C B 0 m=s B lC B B C B @ qg A @ 1:3 kg=m3 ql

L

1000 kg=m3

1 C C C C C; C C C A

1

0

0

ag 0:8 B p C B 105 Pa B C B B C B B ug C B 0 m=s B C ¼B Bu C B 0 m=s B lC B B C B @ qg A @ 1:3 kg=m3 ql

R

1 C C C C C: C C C A

1000 kg=m3

Fig. 3 is obtained at time t = 0.4 s with 100 grids: the results by Zein’s HLLC and present HLLC schemes agree with the exact solution. CFL = 0.9 is used for the computation. The discontinuity of the void fraction is well resolved by the two HLLC schemes but not by the diffusive HLL scheme. Other variables (pressure, velocity, and density) remain constant but not the velocity field shown by the HLL scheme having a small disturbance. In the second problem, the discontinuity of void fraction travels to the right with a constant speed, 10 m/s. The initial left and right conditions are given as follows:

L

1000 kg=m3

1 C C C C C; C C C A

ag 1 0 0:8 B p C B 105 Pa B C B B C B B ug C B 10 m=s B C ¼B Bu C B 10 m=s B lC B B C B @ qg A @ 1:3 kg=m3 0

ql

R

1 C C C C C: C C C A

1000 kg=m3

In Fig. 4, the discontinuity of void fraction is advected to the right for all schemes. Both the Zein’s HLLC and the present HLLC scheme show the same resolution; but the HLL scheme shows poor resolution. No disturbances are observed with other variables that remain constant. These numerical results demonstrate that two HLLC-type schemes satisfy the Abgrall’s principle fairly well. 4.2. Cavitation problem We consider the cavitation shock tube problem simulating the cavitation generation and growth, that has been investigated by several authors [2,6,11]. A 1-m long pipe is initially filled with a two-phase mixture with only 1% dispersed gas at atmospheric pressure. The fluids in the left and in the right part of the midsection are suddenly set into motion in the opposite direction so that cavitation is initiated by the rarefaction waves. The cavitation is

98

G.-S. Yeom, K.-S. Chang / Computers & Fluids 76 (2013) 86–104

Fig. 13. The high-pressure liquid dodecane/low-pressure vapor dodecane shock tube problem with a pressure ratio of 105.

then generated in the center region and grows in time. The initial left and right data are given as follows:

0

ag 1

0

0:01

B p C B 105 Pa B C B B C B B ug C B 100 m=s B C ¼B Bu C B 100 m=s B lC B B C B @ qg A @ 50 kg=m3

ql

L

1000 kg=m3

1

0

ag 1

C C C C C; C C C A

B p C B 105 Pa B C B B C B B ug C B 100 m=s B C ¼B Bu C B 100 m=s B lC B B C B @ qg A @ 50 kg=m3

ql

R

0

0:01

1 C C C C C: C C C A

1000 kg=m3

The EOS parameters used are: cg = 1.4, p1,g = 0 Pa, cl = 1.932, and p1,l = 1.1645  109 Pa. Fig. 5, obtained at time t = 1 ms with 100 grids, shows the numerical results showing the left- and the right-traveling rarefaction waves. The numerical results of the two HLLC schemes are compared with the reference solution that has been obtained earlier by the present authors with HLL scheme on a fine grid of 10,000 cells [6]. In these figures, the present HLLC scheme shows slightly better resolution than the Zein’s HLLC scheme; the HLL scheme with 100 grids performs worst among the three. 4.3. Liquid–gas shock tube problem 4.3.1. High-pressure water/low-pressure air Consider a shock tube filled with water under high pressure at the left, and with air at atmospheric pressure at the right. By sudden rupture of the diaphragm, three waves are launched: a leftgoing rarefaction and two right-going interface and shock waves.

We first solve a water–air shock tube problem with a pressure ratio of 104. A 1-m long tube has a diaphragm at 0.7 m. The initial conditions are as follows:

0

ag 1

0

e

B p C B 10 Pa B C B B C B B ug C B 0 m=s B C ¼B Bu C B 0 m=s B lC B B C B @ qg A @ 50 kg=m3 9

ql

L

1000 kg=m3

1

0

ag 1

C C C C C; C C C A

B p C B 10 Pa B C B B C B B ug C B 0 m=s B C ¼B Bu C B 0 m=s B lC B B C B @ qg A @ 50 kg=m3

0

1e 5

ql

R

1 C C C C C: C C C A

1000 kg=m3

Each fluid is governed by the stiffened-gas EOS, where the parameters used are: cg = 1.4, p1,g = 0 Pa, cl = 4.4, and p1,l = 6  108 Pa. For numerical reasons, a small volume fraction e of the other fluid is presented for in each side of the shock tube to make virtually pure fluids. It can be naturally expected that the lower the value of e, the more accurate the result. In general, the numerical scheme easily becomes unstable as e is lowered to a very small value, typically e = 108. For the present six-equation model, only Paillère et al. [17] solved this problem (a pressure ratio of 104) by the AUSM + scheme, but they used a relatively large value of void fraction, e = 103. In their paper, they said that ‘‘In [25], computations are performed with a valuee = 108. In our case, we have been unable to run the scheme in a satisfactory way with such a small value of void fraction, and have therefore chosen the value e = 103.’’ However, the present HLLC scheme for the sixequation model can solve this problem with a very small value e = 109. Fig. 6, obtained at time t = 229 s with 100 grid cells,

G.-S. Yeom, K.-S. Chang / Computers & Fluids 76 (2013) 86–104

99

Fig. 14. The high-pressure air/low-pressure water shock tube problem with a pressure ratio of 104.

compares the results of the two HLLC schemes with the exact solution. It is noted that the HLL scheme failed to return any solution for this problem. In Fig. 6, both the present and Zein’s HLLC schemes show satisfactory results. Nevertheless the present scheme returns significantly better resolution for the left-traveling rarefaction wave and the right-traveling shock wave. We next investigate the effect of e for the same problem. Fig. 7 shows the numerical results solved by the present HLLC scheme with 5000 grid cells for four different values of e: e = 103, e = 105, e = 107, and e = 109. In this figure, the pressure profiles are plotted on logarithmic scales for clarity. The numerical result solved with e = 103 gives an unacceptable solution, showing a large difference from the exact solution. The results converge to the exact solution as e decreases to e = 109. It is noted that a small difference appears between the results with e = 107 and with e = 109. We now investigate the effect of the EOS parameters. As summarized in Table 1, each of the earlier researchers has been used different parameters of the stiffened-gas EOS for liquid water to solve the water–air shock tube problem. The present HLLC scheme gives robust results with e = 109 for all of the three cases in Table 1 without any difficulties. The results, obtained with 5000 grid cells, are shown in Fig. 8. As shown in this figure, both the wave speed and the wave strength are quite sensitive to the EOS parameters, indicating that an accurate estimation of the EOS parameters for liquid water is very important for this problem.

Finally, we consider a very stiff shock tube problem to show the robustness of the present scheme: initially, the water pressure at the left is set to 1010 Pa (a pressure ratio of 105) instead of 109 Pa, but a volume fraction is kept to the same value e = 109. For the six-equation two-fluid model, such a high pressure ratio has not been solved by the other researchers to the best of our knowledge. Only Saurel et al. [24] solved the problem with an extremely high pressure ratio of 107 using the Baer–Nunziato type model. Note that in [24], they called their model the six-equation model, but it is totally different from the present six-equation model. In Fig. 9, the results, solved by the present HLLC scheme, are plotted at t = 70 s. As grid cells increase, the results converge to the exact solution very well, showing the robustness of the present HLLC scheme.

4.3.2. High-pressure liquid dodecane/low-pressure vapor dodecane In this section, we consider the dodecane liquid–vapor shock tube problem solved by Zein et al. [12]: a 1-m long shock tube is filled with liquid dodecane under high pressure at the left, and with the vapor dodecane at atmospheric pressure at the right. The initial discontinuity is set at 0.75 m. The initial conditions are as follows: Left: pl = 108 Pa, ql = 500 kg/m3, ul = 0 m/s, and ag = e. Right: pg = 105 Pa, qg = 2 kg/m3, ug = 0 m/s, and ag = 1  e.

100

G.-S. Yeom, K.-S. Chang / Computers & Fluids 76 (2013) 86–104

Fig. 15. The high-pressure vapor dodecane/low-pressure liquid dodecane shock tube problem with a pressure ratio of 105.

Fig. 16. Schematic of the solid mixture Hugoniot test problem.

The EOS parameters for liquid and vapor dodecanes are given in Table 2. Fig. 10 compares the present results solved by the modified HLLC scheme with Zein et al.’s results [12] with the same initial data and mesh points as in Zein’s paper. In this figure, the two results are almost the same. Only a small difference appears in the pressure curve at x  0.5 m. This small difference is not due to the numerical scheme but the mathematical model of a two-phase flow: the present paper used the six-equation two-fluid model, while Zein et al. [12] used the Baer–Nunziato type model (the seven-equation model). Next, we apply both the present HLLC scheme and Zein’s HLLC scheme to the present six-equation two-fluid model. The results, obtained with 1000 grid cells, are shown in Fig. 11. In this figure, two schemes shows similar performance, but the present HLLC scheme shows better resolution in the right-traveling shock wave. We now investigate the effect of e for this problem. The present scheme is very robust for a much smaller value e = 108 than

e = 106 used by Zein et al. [12]. In Fig. 12, the results, obtained with 10,000 grid cells, are compared with the exact solution. As can be expected, the numerical results approach to the exact solution more accurately with decreasing e: this is clearly seen in the pressure profile of Fig. 12. In order to show the robustness of the present scheme, we raise pressure of the liquid dodecane at the left to 109 Pa (or 1010 Pa): the resulting pressure ratio is 104 (or 105). The initial discontinuity of a shock tube is set at 0.5 m. The initial conditions are as follows: Left: pl = 109 Pa (or pl = 1010 Pa), ql = 1000 kg/m3, ul = 0 m/s, and

ag = e. Right: pg = 105 Pa, qg = 2 kg/m3, ug = 0 m/s, and ag = 1  e. The numerical results solved by the present HLLC scheme at t = 80 ls are shown in Fig. 13 for pl = 1010 Pa and e = 109. As grid cells increase, the results converge to the exact solution very well for the stiff problem, showing the robustness of the present scheme. Note that in the case of the pressure ratio (pl = 109 Pa to pg = 105 Pa), the numerical results agree well with the exact solution without any difficulty.

4.3.3. High-pressure air/low-pressure water In this section, we consider the reverse configuration: a shock tube filled with air under high pressure at the left, and with water

G.-S. Yeom, K.-S. Chang / Computers & Fluids 76 (2013) 86–104

Fig. 17. Computed pressure profiles solved by the present HLLC scheme for the solid mixture (Cu/Zn) Hugoniot test problem at four different times under piston velocities of 700 m/s and 1000 m/s. (a) u = 700 m/s. (b) u = 1000 m/s.

101

Fig. 18. Comparison of the present results with Andrianov et al. [31] for the solid mixture Hugoniot test problem. Experimental data is compared with computed Hugoniot curves of brass (Cu/Zn). (a) Andrianov et al. [31]. (b) Present result.

at atmospheric pressure at the right. First, we consider a 10 m long shock tube consisting of high-pressure (109 Pa) air at the left and low-pressure (105 Pa) water at the right: a pressure ratio is 104. This problem was solved by Liou et al. [32] by the AUSM+-up scheme. The initial discontinuity is located at 5 m. The initial conditions are as follows:

0

ag

1

BpC B C B C B ug C B C Bu C B lC B C @ Tg A Tl

L

1 1e 9 B 10 Pa C C B C B B 0 m=s C C ¼B B 0 m=s C; C B C B @ 308:15 K A 0

308:15 K

0

ag

1

0

e

1

B 10 Pa C BpC C B B C C B B C B 0 m=s C B ug C C B C ¼B B 0 m=s C: Bu C C B B lC C B B C @ 308:15 K A @ Tg A 5

Tl

R

308:15 K

The EOS parameters used are: cg = 1.4, p1,g = 0 Pa, cl = 2.8, p1,l = 8.5  108 Pa, and Cp,/ = 4186 J/kg K. Liou et al. [32] used e = 107, while we use a smaller value e = 108. It is noted that Zein’s HLLC scheme fails to obtain the solution. Fig. 14 shows the

Fig. 19. The 2D shock–bubble interaction problem.

results solved by the present HLLC scheme, plotted at t = 2.5 ms. As grid cells increase, the results converge to the exact solution very well without any spurious overshoot and oscillation.

4.3.4. High-pressure vapor dodecane/low-pressure liquid dodecane We hereby consider a 10 m long shock tube with the vapor dodecane under very high pressure (1010 Pa) at the left and with liquid dodecane at atmospheric pressure (105 Pa) at the right: a pressure ratio of 105. This stiff problem has not been solved by

102

G.-S. Yeom, K.-S. Chang / Computers & Fluids 76 (2013) 86–104

Fig. 20. Schlieren-type images at the four different times, t = 0.1, 0.3, 0.5, and 0.7 ms. Shock–bubble interaction: (a) by the HLL scheme and (b) by the modified HLLC scheme.

the other researchers. The initial discontinuity of the shock tube is at 3 m. The initial conditions are as follows:

0

ag 1

0

1e

1

B 1010 Pa C BpC C B B C C B B C B 0 m=s C B ug C C B C ¼B B 0 m=s C; Bu C C B B lC C B B C @ 308:15 K A @ Tg A Tl

L

308:15 K

0

ag 1

0

e

1

B 105 Pa C BpC C B B C C B B C B 0 m=s C B ug C C B C ¼B B 0 m=s C: Bu C C B B lC C B B C @ 308:15 K A @ Tg A Tl

R

308:15 K

A small volume fraction is set to a very small value e = 109. Fig. 15 shows the results solved by the present HLLC scheme, plotted at t = 5 ms. It is noted that Zein’s HLLC scheme fails to obtain the solution. As grid cells increase, the results converge to the exact solution very well without any spurious overshoot and oscillation, indicating the robustness of the present HLLC scheme for the stiff shock tube problem. 4.4. Solid mixture Hogoniot test problem In this section, we consider a problem with a very strong shock wave propagating in a mixture of two solid materials, which was investigated by Saurel and Abgrall [11] and later by Andrianov et al. [31] using the Baer–Nunziato type model. This problem has not been, however, solved by the six-equation two-fluid model. Fig. 16 illustrates a schematic of the problem. Brass, a mixture of copper and zinc, is considered as a base material. Both copper and zinc are treated as compressible medium governed by the stiffened-gas EOS; the parameters are: ccopper = 4.22, p1,copper = 32.32  109 Pa, czinc = 4.17, and p1,zinc = 15.71  109 Pa. The initial conditions are as follows:

faster as the piston velocity increases. For several piston speeds, the corresponding shock speeds are calculated by the present HLLC scheme using the present six-equation model. Fig. 18 shows the comparison of the present results with the results by Andrianov et al. [31], where the experimental data is given in [33]. The present result agrees well with both the experimental data and the results by Andrianov et al. 4.5. Two-dimensional shock–bubble interaction problem Finally we solve the two-dimensional shock–bubble interaction problem considered earlier [27,28]. This problem treats collision of a shock wave with a circular R22 gas bubble in air. Both the air and the R22 gas are assumed governed by the stiffenedgas EOS; the parameters are: cair = 1.4, p1,air = 0 Pa, cR22 = 1.249, and p1,R22 = 0 Pa. Fig. 19 shows the problem geometry: the shock tube is 350 mm long and 89 mm wide, the circular R22 gas bubble has a radius of 25 mm, and the shock wave initially located at 5 mm from the right boundary has a shock Mach number of 1.22. The initial pre- and post-shock condition in air are given by

0 1 p BuC B C B C @v A

q

0

1:59  105 Pa

B 113:5 m=s B ¼B @ 0 m=s

1 C C C; A

1:686 kg=m3 1 1:01325  105 Pa C B 0 m=s C B ¼B C: A @ 0 m=s 0

preshock

0 1 p BuC B C B C @v A

q

postshock

1:225 kg=m3 The initial condition inside the R22 gas bubble is

Copper: pcopper = 105 Pa, qcopper = 8924 kg/m3, ucopper = 0 m/s, and acopper = 0.71. Zinc: pzinc = 105 Pa, qzinc = 7139 kg/m3, uzinc = 0 m/s, and azinc = 0.29. A piston boundary condition on the left wall is imposed to initiate the shock wave. Fig. 17 shows the pressure profiles for two piston speeds of 700 m/s and 1000 m/s at four different times. It is clearly observed that the shock wave becomes stronger and

0 1 p BuC B C B C @v A

q

0 B B ¼B @ R22

1 1:01325  105 Pa C 0 m=s C C: A 0 m=s 3:863 kg=m3

In order to represent virtually pure fluids using the two-fluid model, a small amount (0.1%) of air and R22 gas are added respectively to the R22 gas bubble and to the air.

G.-S. Yeom, K.-S. Chang / Computers & Fluids 76 (2013) 86–104

103

Fig. 21. 3D plots of the void fraction for the shock–bubble interaction, given at the four different times, t = 0.1, 0.3, 0.5, and 0.7 ms. Left column is the HLL scheme and right column is the modified HLLC scheme.

Fig. 20 compares, at the four different time instants, t = 0.1, 0.3, 0.5, and 0.7 ms, the Schileren-type image of the air density solved by the earlier HLL scheme [7] with that solved by the present HLLC scheme. We extended both HLL and HLLC schemes to the secondorder for this computation. We used 779,744 grid cells and CFL = 0.4. In this figure, the present HLLC scheme presents fine resolution relative to the HLL scheme with respect to the material interface, transmissive and reflected shock waves. Fig. 21 shows the 3D material interface between the air and the R22 gas bubble represented by the void fraction at the same four different time instants. The present HLLC scheme returns very sharp material interface in contrast to the quite smeared results of the HLL scheme.

5. Conclusions We have developed a modified HLLC-type Riemann solver for the compressible six-equation two-fluid model. The modified HLLC scheme restores the missing characteristic fields of the recent HLLC scheme by Zein et al. [12]. The modified HLLC scheme has been tested on a few one-dimensional two-phase shock tube problems and on a two-dimensional shock–bubble interaction problem. The numerical results indicate that the present HLLC scheme applied to the six-equation model is, in fact, more accurate and satisfactory than the recent HLLC scheme: these are summarized in Table 3. It should be noted that the both present and Zein’s HLLC schemes work well for the two-phase flows where the drag force

104

G.-S. Yeom, K.-S. Chang / Computers & Fluids 76 (2013) 86–104

Table 3 Summary of the present (bold face) and earlier results for the liquid–gas shock tube problem. Model

Authors

Scheme

Fluids (HP/LP)a

PRb

e

Baer– Nunziato type

Saurel and Abgrall [11] Andrianov et al. [31] Saurel et al. [29] Saurel et al. [24] Zein et al. [12]

HLL

Water/air

104

108

VFRoe

Water/air

103

108

Rel-Projc

Dodecane liquid/vapor Water/air

103

108

107

106

HLLC

Dodecane liquid/vapor

103

106

Paillère et al. [17] Chang and Liou [30]

AUSM+

Water/air

104

103

AUSM+up

Water/air

103

105

Liou et al. [32]

AUSM+up new AUSM+up Modified HLLC

Air/water Air/water

104 104

107 107

Water/air

2

106

Water/air

105

109

Dodecane liquid/vapor Air/water Dodecane vapor/liquid

105

109

104 105

108 109

Six-equation

Niu et al. [5]

Present

a b c

HLLC

HP: high-pressure; LP: low-pressure. Pressure ratio. Relaxation-projection method.

between phases is very large so that the two phasic velocities become equal very quickly. Otherwise, the both HLLC-type schemes may produce numerical oscillation and fail to converge. The future work will be the extension of the present HLLC scheme for the twophase flows with a small drag force. References [1] Ishii M, Hibiki T. Thermo-fluid dynamics of two-phase flow. New York: Springer; 2006. [2] Saurel R, Lemetayer O. A multiphase model for compressible flows with interfaces, shocks, detonation waves and cavitation. J Fluid Mech 2001;431:239–71. [3] Toumi I. An upwind numerical method for two-fluid two-phase flow models. Nucl Sci Eng 1996;123:147–68. [4] Städtke H, Franchello G, Worth B. Numerical simulation of multi-dimensional two-phase flow based on flux vector splitting. Nucl Eng Des 1997;177:199–213. [5] Niu YY, Lin YC, Chang CH. A further work on multi-phase two-fluid approach for compressible multi-phase flows. Int J Numer Meth Fluids 2008;58:879–96. [6] Yeom GS, Chang KS. numerical simulation of two-fluid two-phase flows by HLL scheme using an approximate Jacobian matrix. Numer Heat Transfer B-Fund 2006;49:155–77. [7] Yeom GS, Chang KS. Two-dimensional two-fluid two-phase flow simulation using an approximate Jacobian matrix for HLL scheme. Numer Heat Tr B-Fund 2009;56:372–92.

[8] Yeom GS, Chang KS. Wave characteristics of two-phase flows predicted by HLL scheme using interfacial friction terms. Numer Heat Tr A-Appl 2010;58:356–84. [9] Yeom GS, Chang KS. Shock wave diffraction about a wedge in a gas– microdroplet mixture. Int J Heat Mass Transfer 2010;53:5073–88. [10] Baer MR, Nunziato JW. A two-phase mixture theory for the deflagration-todetonation transition (DDT) in reactive granular materials. Int J Multiphase Flow 1986;12:861–89. [11] Saurel R, Abgrall R. A multiphase Godunov method for compressible multifluid and multiphase flows. J Comput Phys 1999;150:425–67. [12] Zein A, Hantke M, Warnecke G. Modeling phase transition for compressible two-phase flows applied to metastable liquids. J Comput Phys 2010;229:2964–98. [13] Toro EF. Riemann solvers and numerical methods for fluid dynamics: a practical introduction. 2nd ed. Berlin: Springer-Verlag; 1999. [14] Sainsaulieu L. Finite volume approximation of two phase-fluid flows based on an approximate roe-type riemann solver. J Comput Phys 1995;121:1–28. [15] Cortes J. On the construction of upwind schemes for non-equilibrium transient two-phase flows. Comput Fluids 2002;31:159–82. [16] Coquel F, El Amine K, Godlewski E, Perthame B, Rascle P. A numerical method using upwind schemes for the resolution of two-phase flows. J Comput Phys 1997;136:272–88. [17] Paillère H, Corre C, Garcia Cascales JR. On the extension of the AUSM+ scheme to compressible two-fluid models. Comput Fluids 2003;32:891–916. [18] Harten A, Lax PD, van Leer B. On upstream differencing and Godunov-type schemes for hyperbolic conservation laws. SIAM Rev 1983;25:35–61. [19] Toro EF, Spruce M, Speares W. Restoration of the contact surface in the HLLRiemann solver. Shock Waves 1994;4:25–34. [20] Tokareva SA, Toro EF. HLLC-type Riemann solver for the Baer–Nunziato equations of compressible two-phase flow. J Comput Phys 2010;229:3573–604. [21] Barre F, Bernard M. The CATHARE code strategy and assessment. Nucl Eng Des 1990;124:257–84. [22] Tiselj I, Petelin S. Modelling of two-phase flow with second-order accurate scheme. J Comput Phys 1997;136:503–21. [23] Deledicque V, Papalexandris MV. An exact Riemann solver for compressible two-phase flow models containing non-conservative products. J Comput Phys 2007;222:217–45. [24] Saurel R, Petitpas F, Berry RA. Simple and efficient relaxation methods for interfaces separating compressible fluids, cavitating flows and shocks in multiphase mixtures. J Comput Phys 2009;228:1678–712. [25] Toro EF. A weighted average flux method for hyperbolic conservation laws. Proc Roy Soc A 1989;423:401–18. [26] Abgrall R. How to prevent pressure oscillations in multicomponent flow calculations: a quasi conservative approach. J Comput Phys 1996;125:150–60. [27] Quirk JJ, Karni S. On the dynamics of a shock–bubble interaction. J Fluid Mech 1996;318:129–63. [28] Shyue KM. A wave-propagation based volume tracking method for compressible multicomponent flow in two space dimensions. J Comput Phys 2006;215:219–44. [29] Saurel R, Petitpas F, Abgrall R. Modelling phase transition in metastable liquids: application to cavitating and flashing flows. J Fluid Mech 2008;607:313–50. [30] Chang CH, Liou MS. A robust and accurate approach to computing compressible multiphase flow: stratified flow model and AUSM+-up scheme. J Comput Phys 2007;225:840–73. [31] Andrianov N, Saurel R, Warnecke G. A simple method for compressible multiphase mixtures and interfaces. Int J Numer Meth Fluids 2003;41:109–31. [32] Liou MS, Chang CH, Nguyen L, Theofanous TG. How to solve compressible multifluid equations: a simple, robust, and accurate method. AIAA J 2008;46:2345–56. [33] Marsh SP. LASL shock Hugoniot data. Berkely: University of California Press; 1980. [34] Niu YY. Advection upwinding splitting method to solve a compressible twofluid model. Int J Numer Meth Fluids 2001;36:351–71.