A gas kinetic scheme for the Baer–Nunziato two-phase flow model

Journal of Computational Physics 231 (2012) 7518–7536

Contents lists available at SciVerse ScienceDirect

Journal of Computational Physics journal homepage: www.elsevier.com/locate/jcp

A gas kinetic scheme for the Baer–Nunziato two-phase flow model q Liang Pan a,c, Guiping Zhao d, Baolin Tian a,b, Shuanghu Wang a,b,⇑ a

Institute of Applied Physics and Computational Mathematics, Beijing 100088, China Key Laboratory of Computational Physics, Institute of Applied Physics and Computational Mathematics, Beijing 100088, China The Graduate School of China Academy of Engineering Physics, Beijing 100088, China d National Natural Science Foundation of China, Beijing 100085, China b c

a r t i c l e

i n f o

Article history: Received 8 September 2011 Received in revised form 11 April 2012 Accepted 18 April 2012 Available online 31 May 2012 Keywords: Baer–Nunziato model Generalized BGK model Gas kinetic method

a b s t r a c t Numerical methods for the Baer–Nunziato (BN) two-phase flow model have attracted much attention in recent years. In this paper, we present a new gas kinetic scheme for the BN two-phase flow model containing non-conservative terms in the framework of finite volume method. In the view of microscopic aspect, a generalized Bhatnagar– Gross–Krook (BGK) model which matches with the BN model is constructed. Based on the integral solution of the generalized BGK model, we construct the distribution functions at the cell interface. Then numerical fluxes can be obtained by taking moments of the distribution functions, and non-conservative terms are explicitly introduced into the construction of numerical fluxes. In this method, not only the complex iterative process of exact solutions is avoided, but also the non-conservative terms included in the equation can be handled well. Ó 2012 Elsevier Inc. All rights reserved.

1. Introduction The compressible two-phase flows arise in various natural phenomena and technological applications. One popular approach to describe two-phase flows is to use the so-called homogenized or averaged mixture models. In the framework of these models, the mixture is assumed to consist of two separate, identifiable and interpenetrating continua that are in thermodynamic non-equilibrium with each other. Examples of mixture models for compressible two-phase flows are the Baer–Nunziato (BN) model of Baer and Nunziato [3], as well as those of Saurel and Abgrall [18], Bdzil et al. [4] and Papalexandris [17], etc. Most of the two-phase flow models contain non-conservative products which are part of the expressions that describe momentum and energy transfer between two phases. The presence of the non-conservative terms introduces serious analytical and computational difficulties. It is the non-conservative terms that do not permit the governing equations to be written in divergence form. Concerning the mathematical analysis, it is clear that the solution may become discontinuous, so that the weak solution needs to be introduced. However, since the governing equations is non-conservative, one cannot use the corresponding definition from the theory of distributions which is used for conservation laws. As a result, the Rankine– Hugoniot shock conditions cannot be defined as it is done for conservation laws. There has been recent progresses in the mathematical analysis of non-conservative hyperbolic equations Dal Maso et al. [16] and Lefloch and Tzavaras [12].

q

This work was supported by the Natural Science Foundation of China (NSFC) Nos. 10931004, 11171037 and 91130021.

⇑ Corresponding author at: Institute of Applied Physics and Computational Mathematics, Beijing 100088, China.

E-mail addresses: [email protected] (L. Pan), [email protected] (G. Zhao), [email protected] (B. Tian), [email protected] (S. Wang). 0021-9991/$ - see front matter Ó 2012 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.jcp.2012.04.049

L. Pan et al. / Journal of Computational Physics 231 (2012) 7518–7536

7519

Because of these serious analytical difficulties, the derivation of efficient algorithms for the numerical treatment of the governing equations becomes a very challenging task and one has difficulties in dealing with the non-conservative terms. Additionally, since the non-conservative products do not represent fluxes, they have to be integrated as source terms. Currently, the theory of numerical methods for non-conservative systems is largely absent. There is no criterion how one should design a numerical scheme for such systems. Recently, there has also been progress in the construction of numerical methods for the flow models of interest. For example, Saurel and Abgrall [18], Saurel and Lemetayer [13] and Gavrilyuk and Saurel [10] presented numerical approximations of the nozzling terms based on a condition of Abgrall [1]. Since its introduction, the Baer–Nunziato model [3] for two-phase flows has received much attention. A mathematical analysis of the structure of the governing equations, including its Riemann invariants and simple wave solutions and the description of the wave fields of the hyperbolic system, was described by Embid and Baer [9]. By making use of the jump conditions across the solid contact, Andrianov and Warnecke [2] constructed inverse exact solutions corresponding to prescribed states on either side of the solid contact. Schwendeman et al. [19] provided the exact Riemann solver of the flow models and a numerical procedure. They described a two-stage iterative procedure to obtain exact solutions of the Riemann problem directly for left and right states. Deledicque and Papalexandris [8] proposed a different procedure for obtaining numerical solutions of the Riemann problem. The algorithm proposed is based on the possible solutions into four principal configurations. The admissibility of each configuration is examined sequentially, and the intermediate states are calculated via two secant iterations. Approximate Riemann solvers are also presented [21]. In recent years, the gas kinetic scheme based on the Bhatnagar–Gross–Krook (BGK) model [5,7,6] for compressible fluids proposed in [23–25,30] has attracted much attention. Based on the Chapmann-Enskog expansion, Euler equation as well as Navier–Stokes equation can be derived from the BGK model. In the framework of finite volume method, the BGK scheme makes use of the local integral solution of the collisional BGK model at a cell interface to compute a time-dependent gas distribution function from which the numerical fluxes can be obtained in the gas evolution stage. As the BGK model is a statistical model, particle transport and collision are coupled in the whole gas evolution process, and the particle collision time controls physical dissipative coefficients in the macroscopic equations. Since the gas evolution is associated with a relaxation process i.e. from a non-equilibrium state to an equilibrium one, the entropy condition is satisfied automatically by the BGK scheme. Once the physical structure can be well resolved by the numerical cell size, in smooth regions, the scheme automatically gives an accurate compressible Navier–Stokes solution; in the discontinuous regions, the BGK scheme generates a stable and crisp shock transition because of the delicate dissipative mechanism. The BGK scheme has also been extended to multicomponent flow [11,15,26], magnetohydrodynamics [28,20], hyperbolic conservation laws with source terms [23] and rarefired gas flow [29]. And the high order BGK schemes have also been also developed [27,14]. Recently, She and Zhao [32] proposed a gas kinetic method for the BN two-phase flow model containing non-conservative terms. In view of microscopic aspect, a collisionless gas kinetic model which matches with the BN model is constructed, then a Kinetic Flux Vector Splitting (KFVS) scheme is proposed. Through this method, the complex iterative process [2,19,8] can be avoided. However in the gas evolution stage, KFVS scheme is based on the collisionless model [23]. The collisions of the particles are not taken into consideration, so a physical model in which particle collisions in the gas evolution stage are considered should be constructed. In [31], a generalized BGK model with the inclusion of the gravitational force was proposed, and the gas-kinetic BGK scheme was extended to the shallow-water equations with source terms, where the effect of the gravitational force was taken into account. In this paper, we make use of the generalized BGK model [31] in which the non-conservative term is considered as the external forces for the BN model, and we take the particle collisions into account which are not considered in the collisionless gas kinetic model. Similar with the conventional gas kinetic methods, the distribution functions, which is based on the integral solution of the generalized BGK model, is constructed at the cell interface. Then numerical fluxes can be got by taking moments of the distribution functions, and non-conservative terms are introduced into the evolution and construction of numerical fluxes. Compared with the exact and approximate Riemann solves, in our scheme, we avoid the iterative process and the non-conservative term included in the equation can be handled well. The collisions of the particles in the evolution stage are taken into account. This paper is organized as follows. In Section 2, we will introduce the Baer–Nunziato model and the generalized BGK model, and show the relationship between them. In Section 3, the outline of the numerical scheme will be given. In Section 4, some numerical examples are presented to validate the scheme and we will give some discussions about the current scheme. The last section is the conclusion. 2. The Baer–Nunziato model and the generalized BGK model 2.1. The Baer–Nunziato model The Baer–Nunziato (BN) model under consideration treats each phase as two separate and co-existing continua. For the sake of simplicity in notation we will refer to the carrier phase as the gaseous phase and to the particulate one as the solid phase. In this paper, we only consider the homogeneous Baer–Nunziato model. The governing equations can be written in the following form:

7520

L. Pan et al. / Journal of Computational Physics 231 (2012) 7518–7536

@/s @/ þ us s ¼ 0; @t @x @/s qs @/s qs us þ ¼ 0; @t @x @/s qs us @/s ðqs u2s þ ps Þ @/ ¼ pg s ; þ @x @t @x @/s qs Es @/s us ðqs Es þ ps Þ @/ ¼ pg us s ; þ @x @t @x @/g qg @/g qg ug þ ¼ 0; @t @x @/g qg ug @/g ðqg u2g þ pg Þ @/ þ ¼ pg s ; @t @x @x @/g qg Eg @/g ug ðqg Eg þ pg Þ @/ ¼ pg us s ; þ @x @t @x

ð1:aÞ ð1:bÞ ð1:cÞ ð1:dÞ ð1:eÞ ð1:fÞ ð1:gÞ

where each phase is assigned a density qa , pressure pa , specific internal energy ea , velocity ua , total specific energy Ea ¼ ea þ 12 u2a and volume fraction /a , where a ¼ g; s, and g; s denote the carrier gaseous phase and the particulate solid phase, respectively. The volume fractions represent the percentages of volume occupied by the constituents. It is also assumed that the following saturation condition holds:

/g þ /s ¼ 1: The above system is closed with the equation of state (EOS) for each phase. In order to avoid difficulties with thermodynamic modeling, both the ideal gas EOS:

ea ¼

pa

qa ðca  1Þ

;

ð2Þ

and the stiffened gas EOS:

ea ¼

pa þ pa ; qa ðca  1Þ

ð3Þ

could be used, where ca are the ratio of the specific heats, pa are constants, a ¼ s; g. To simplify the governing equations, we denote qð1Þ ¼ /s qs ; qð2Þ ¼ /g qg ; Pð1Þ ¼ /s ps , Pð2Þ ¼ /g pg ; Eð1Þ ¼ Es ; Eð2Þ ¼ Eg ; U ð1Þ ¼ us , U ð2Þ ¼ ug , where the solid and gaseous phase are represented by 1 and 2, respectively. Eq. (1.a) can be written into the following conservative form

@ qð1Þ =/s @ qð1Þ U ð1Þ =/s þ ¼ 0: @t @x

ð4Þ

Let Z ¼ 1=/s , we have the following governing equations:

@ qð1Þ Z @ qð1Þ U ð1Þ Z þ ¼ 0; @t @x ð1Þ @ qð1Þ @ qð1Þ U þ ¼ 0; @t @x ð1Þ ð1Þ ð1Þ @q U @ðq ðU ð1Þ Þ2 þ Pð1Þ Þ @/ ¼ pg s ; þ @x @t @x @ qð1Þ Eð1Þ @U ð1Þ ðqð1Þ Eð1Þ þ Pð1Þ Þ @/ ð1Þ ¼ U pg s ; þ @x @t @x @ qð2Þ @ qð2Þ U ð2Þ þ ¼ 0; @t @x @ qð2Þ U ð2Þ @ðqð2Þ ðU ð2Þ Þ2 þ Pð2Þ Þ @/ ¼ pg s ; þ @x @t @x @ qð2Þ Eð2Þ @U ð2Þ ðqð2Þ Eð2Þ þ Pð2Þ Þ @/ ¼ U ð1Þ pg s : þ @x @t @x

ð5:aÞ ð5:bÞ ð5:cÞ ð5:dÞ ð5:eÞ ð5:fÞ ð5:gÞ

We can observe that the system Eqs. (1) and (5) cannot be written in divergence form, i.e., it is non-conservative. Indeed, the non-conservative terms prevent it from being written in divergence form. Note that for the case /a ¼ const, the system decouples into the two sets of Euler equations for each phase. For the mixture, the conservative equations can be obtained by summing the corresponding single-phase equations.

L. Pan et al. / Journal of Computational Physics 231 (2012) 7518–7536

7521

2.2. The generalized BGK model In this section, we propose a generalized BGK model to describe the BN two-phase flow model in one dimensional case. In this model, the non-conservative terms describe the interaction between the solid and gaseous phases, so that the non-conservative terms can be treated as the external force in the generalized BGK model for each phase. The distribution functions are functions of spatial coordinate x, velocity-space coordinate u and time t. The generalized BGK model with inclusion of external force may be written as: ð1Þ

ft

ð2Þ

ft

ð1Þ

þ uf x þ xfuð1Þ ¼ ð2Þ

þ uf x  xfuð1Þ ¼

g ð1Þ  f ð1Þ

s g ð2Þ  f ð2Þ

s

ð6Þ

;

ð7Þ

; ðiÞ

where the solid phase and the gaseous phase are represented by 1 and 2, respectively. f is the distribution function of the solid or gaseous phases and g ðiÞ is the equilibrium state approaching f ðiÞ . s is the particle collision time to be specified and x is considered as the external force. The equilibrium states g ð1Þ ; g ð2Þ are all local Maxwellians in velocity space, which can be expressed as follows:

g ð1Þ ðx; v ; tÞ ¼ qð1Þ g ð2Þ ðx; v ; tÞ ¼ qð2Þ

 K 12þ1 k1

p  K 22þ1 k2

p

ek1 ððv U

ð1Þ 2

ek2 ððv U

ð2Þ 2

Þ þn21 Þ

Þ þn22 Þ

; ;

mi where ki is function of temperature, Boltzmann’s constant and molecule mass, with the relationship ki ¼ 2kT . i   K are the component of the internal particle velocity in Ki dimensions, and ni ¼ n1i ; n2i ; . . . ; ni i  2  2  2 K n2i ¼ n1i þ n2i þ . . . þ ni i . The ratio of the principal specific heats denoted by ci . In 1-D case , we can obtain the follow-

ci ing relationship K i ¼ 3 c 1. U i is the corresponding macroscopic flow velocity. Due to the mass, momentum and energy coni

servation, the collision term satisfies the compatibility condition

Z

g ðiÞ  f ðiÞ

s

WðiÞ a du dni ¼ 0;

ð8Þ K

2 1 1 2 i where WðiÞ a ¼ ð1; u; 2 ðu þ ni ÞÞ, where a ¼ 1; 2; 3, du dni ¼ du dni . . . dni .

2.3. Relationship between two models Derivations of the Euler and Navier–Stokes Equations from the Boltzmann equation can be found in Chapmann and Cowling [7] and from the BGK equation can be found in Cercignani [6] and Xu [23]. In this section, we present the governing equations of the BN two-phase model derived from the generalized BGK model Eqs. (6) and (7). Denote

  Sð1Þ ¼ x 0; qð1Þ ; U ð1Þ qð1Þ ;

ð9Þ

where

x¼

pg @/s Pð1Þ @/s ¼ : ð1Þ q @x ð1  /s Þqð1Þ @x

ð10Þ

For a local equilibrium state f ð1Þ ¼ g ð1Þ and f ð2Þ ¼ g ð2Þ , the governing equations of the BN model can be obtained by taking the 2 1 2 moment of WðiÞ a ¼ ð1; u; 2 ðu þ ni ÞÞ to the generalized BGK model Eqs. (6) and (7). This yields

Z

and

Z

ð1Þ

ð11Þ

ð2Þ

ð12Þ

ð1Þ ðg t þ ug xð1Þ þ xg ð1Þ u ÞWa du dn1 ¼ 0

ð2Þ ðg t þ ug xð2Þ  xg ð1Þ u ÞWa du dn2 ¼ 0:

Consider the relationship between macroscopic variables and the local equilibrium state, we have:

W ðiÞ ¼

Z

FðW ðiÞ Þ ¼

ðiÞ WðiÞ a g du dni ;

Z

ðiÞ uWðiÞ a g du dni ;

7522

L. Pan et al. / Journal of Computational Physics 231 (2012) 7518–7536

SðiÞ ¼ ð1Þi

Z

ðiÞ uWðiÞ a g u du dni ;

  where W ¼ ðq; qU; qEÞ; FðWÞ ¼ qU; qU 2 þ P; UðqE þ PÞ . With the compatibility condition Eq. (8), we can obtain the macroscopic governing equations Eqs. (5.b)–(5.d): ð1Þ W ð1Þ Þx ¼ Sð1Þ : t þ FðW

R ð1Þ In Eq. (12), g ð1Þ contains n21 term which is different from the integral over du dn2 . In the following, xg u Wð2Þ a du dn2 will be R ð1Þ ð1Þ approximately considered as xg u Wa du dn1 , then we could obtain the macroscopic governing equations (5.e)–(5.g): ð2Þ W ð2Þ Þx ¼ Sð1Þ : t þ FðW

Multiplying Eq. (6) by Z ¼ /1s and taking the moment of W1 ¼ 1, we can get Eq. (5.a). Thus the relationship between the BN two-phase model and the generalized BGK model has been presented.It could be observed that the governing equations of the BN model is only the zeroth order result of Chapmann–Enskog expansion. 3. Numerical method 3.1. The finite volume gas kinetic scheme We have discussed the relationship between the kinetic theory and macroscopic equations. In this section, the numerical scheme will be presented. A uniform mesh with xj ¼ jh; ðj ¼ 0; 1; 2; . . .Þ is considered and h is the mesh size. In order to de2 1 2 velop a finite volume gas-kinetic scheme, taking moments of WðiÞ a ¼ ð1; u; 2 ðu þ ni ÞÞ in Eqs. (6) and (7),integrate with respect n nþ1 to du dn in phase space, dx dt in X ¼ ½xj1=2 ; xjþ1=2   ½t ; t , we have:

Z Z Z

ð1Þ

ðft

X

ð1Þ

þ uf x þ xfuð1Þ ÞWð1Þ a du dn1 dx dt ¼

Z Z Z

g ð1Þ  f ð1Þ

s

X

Wð1Þ a du dn1 dxdt;

with the compatibility condition Eq. (8), we obtain: nþ1;ð1Þ

Wj and

Z Z Z

n;ð1Þ

¼ Wj

ð2Þ

ðft

X

þ

1 Dx

Z

t nþ1

tn

ð1Þ

n;ð1Þ

ð1Þ

ðF j1=2 ðtÞ  F jþ1=2 ðtÞÞdt þ DtSj

ð2Þ

þ uf x  xfuð1Þ ÞWð2Þ a du dn2 dxdt ¼

Z Z Z

ð13Þ

;

g ð2Þ  f ð2Þ

s

X

Wð2Þ a du dn2 dx dt:

Also we have: nþ1;ð2Þ

Wj

n;ð2Þ

¼ Wj

þ

1 Dx

Z

t nþ1

tn

ð2Þ

ð2Þ

n;ð1Þ

ðF j1=2 ðtÞ  F jþ1=2 ðtÞÞdt  DtSj

ð14Þ

;

n;ðiÞ

ðiÞ

is the average value of W ðiÞ in the j-th cell at t ¼ t n ; F jþ1=2 ðtÞ is the numerical flux across the cell inter  @/ n;ð1Þ Pn;ð1Þ s ¼ xn;ð1Þ 0; qn;ð1Þ ; U n;ð1Þ qn;ð1Þ and xn;ð1Þ ¼ . n;ð1Þ n;ð1Þ @x

where Dx ¼ h and W j n;ð1Þ

face x ¼ xjþ1=2 , Sj

ð1/s

Þq

3.2. The BGK fluxes In the finite-volume gas-kinetic scheme, the numerical fluxes are needed at the cell interfaces. With the relationship between distribution function and macroscopic variables and numerical fluxes, the time-dependent gas distribution functions need to be constructed at the cell interface. Before we construct the gas-kinetic scheme, let us consider the particle velocity change due to the external force and its modification to the gas distribution function. Suppose that a particle is moving from to x0 to x in a short time interval t  t0 . Because of the acceleration, the particle velocity will change from to u0 to u ¼ u0  xðt  t0 Þ. So due to the particle velocity change, an equilibrium state g 0 at time t 0 can be expressed as:

 Kþ1  Kþ1 k 2 k½ðu0 UÞ2 þn2  k 2 k½ðuUÞ2 þn2  k½2ðuUÞxðtt0 Þþxðtt0 Þ2  0 0 2 g 0 ðx0 ; t0 ; u0 ; nÞ ¼ q e ¼q e e ¼ g  ek½2ðuUÞxðtt Þþxðtt Þ  ;

p

p

where g is the equilibrium state at time t. Using the Taylor expansion and neglect the ðt  t0 Þ2 order term, the above equilibrium state can be expressed as:

g 0 ðx0 ; t0 ; u0 ; nÞ ’ g½1 þ 2xðu  UÞðt  t0 Þ:

ð15Þ

With the cell averaged mass, momentum, energy of each phase i ¼ 1; 2 in the j-th cell denoted by ðiÞ ðiÞ ðiÞ ðiÞ ðiÞ ðiÞ W j ¼ ðqj ; qj U j ; qj Ej Þ, the reconstruction is needed to prepare the initial data. For the gas kinetic scheme, the reconstruc-

7523

L. Pan et al. / Journal of Computational Physics 231 (2012) 7518–7536

tion techniques are directly applied to the conservative variables. To second order accuracy, the interpolated value in the j-th ðiÞ cell is denoted by W j ðxÞ. The interpolated value in the j-th cell, xj1=2 < x < xjþ1=2 , can be written as:

  ðiÞ ðiÞ ðiÞ ðiÞ W j ðxÞ ¼ W j þ L sjþ ; sj ðx  xj Þ;

where Lðsjþ ; sj Þ is a nonlinear limiter, with sjþ ¼ ðW jþ1  W j Þ=h and sj ¼ ðW j  W j1 Þ=h. The conservative variables can be got at two end points W ðiÞ ðxj1=2 Þ; W ðiÞ ðxjþ1=2 Þ in the j-th cell after the reconstruction. For example, the MUSCL limiter is:

  1 Lðsjþ ; sj Þ ¼ Sgnðsjþ ; sj Þmin jsjþ þ sj j; jsjþ j; jsj j ; 2 and the Van Leer limiter is:

2jsjþ jjsj j ; jsjþ j þ jsj j

Lðsjþ ; sj Þ ¼ Sgnðsjþ ; sj Þ

where Sgnð; Þ is the sign function. In the following, we will construct the numerical fluxes at the cell interfaces for each phase. For the solid phase, by solving the BGK model Eq. (6), the integral solution can be written in the following form:

f ð1Þ ðxjþ1=2 ; t; u; n1 Þ ¼

1

Z

s

t

0

0

0

ð1Þ

g ð1Þ ðx0 ; t 0 ; u0 ; n1 Þeðtt Þ=s dt þ et=s f0 ðxjþ1=2  utÞ;

ð16Þ

ð1Þ

where f0 is the real gas distribution function of f ð1Þ at t ¼ 0 and g ð1Þ is the corresponding equilibrium state. xjþ1=2 is the location of the interface and xjþ1=2 ¼ x0 þ u0 ðt  t0 Þ þ 12 xðt  t 0 Þ2 and u ¼ u0 þ xðt  t 0 Þ. ð1Þ Based on the reconstructed data, we can construct the initial gas distribution function f0 and the equilibrium state g ð1Þ around the cell interface:

  ð1Þ ð1Þ  ð1Þ  f0 ¼ g l0 1 þ a1l ðx  xjþ1=2 Þ ð1  Hðx  xjþ1=2 ÞÞ þ g r0 1 þ a1r ðx  xjþ1=2 Þ Hðx  xjþ1=2 Þ

ð17Þ

ð1Þ

g ð1Þ ¼ g 0 ð1 þ ð1  Hðx  xjþ1=2 ÞÞa1r ðx  xjþ1=2 Þ þ Hðx  xjþ1=2 Þa1l ðx  xjþ1=2 ÞÞ ð1Þ ð1Þ g l0 ; g r0

where interface:

ð1Þ g l0

ð1Þ g r0

ð1Þ g0

and

ð1Þ l

¼q

ð1Þ r

¼q

ð1Þ 0

¼q

ð1Þ g0

k1l

are local Maxwellian distribution functions located to the left, to the right and in the middle of a cell

!K 12þ1

p

1 2

þn21 Þ

;

1

1 2

þn21 Þ

;

1

1 2

þn21 Þ

:

!K 12þ1

ekr ððuUr Þ

p k10

1

ekl ððuUl Þ

p k1r

ð18Þ

!K 12þ1

ek0 ððuU0 Þ

The dependence of the corresponding slopes a1l ; a1r , a1r ; a1l on the particle velocities is obtained from the Taylor expansion of a Maxwellian and have the following form:

 1 2 u þ n21 ; 2  1 a1r ¼ a11r þ a12r u þ a13r u2 þ n21 ; 2  1 a1l ¼ a11l þ a12l u þ a13l u2 þ n21 ; 2  1 a1r ¼ a11r þ a12r u þ a13r u2 þ n21 ; 2 a1l ¼ a11l þ a12l u þ a13l

ð1Þ

ð1Þ

where all the corresponding coefficients are local constants. With the above definition, the equilibrium states g l0 and g r0 can be determined from macroscopic point-wise values at the cell interface:

Z

ð1Þ ð1Þ Wð1Þ a g l0 du dn1 ¼ W j ðxjþ1=2 Þ; ð1Þ

Z

ð1Þ ð1Þ Wð1Þ a g r0 du dn1 ¼ W jþ1 ðxjþ1=2 Þ:

The parameters a1l and a1r in f0 can be obtained by using relationship between the distribution function and the slopes of macroscopic variables on the left and right hand sides:

7524

L. Pan et al. / Journal of Computational Physics 231 (2012) 7518–7536

Z Z

ð1Þ

1 ð1Þ Wð1Þ a al g l0 du dn1 ¼

1 ð1Þ Wð1Þ a ar g r0 du dn1 ¼

ð1Þ

2ðW j ðxjþ1=2 Þ  W j Þ h ð1Þ 2ðW jþ1



h

The slope a1l can be computed from: ð1Þ

0

ð1Þ

2ðW j ðxjþ1=2 Þ  W j Þ

¼

qð1Þ l h

;

ð1Þ W jþ1 ðxjþ1=2 ÞÞ

a11l B 1 M 1l ab @ a2l a13l

:

1 C A:

R ð1Þ ð1Þ ð1Þ 1 where the matrix M 1l Wa Wb g l0 du dn1 and the solution of ða11l ; a12l ; a13l Þ can be found in [23,24]. Similarly, the slope a1r ab ¼ qð1Þ l can be obtained. Eq. (17) only shows the gas distribution function around the interface xjþ1=2 at t ¼ 0, we need to figure out the particle arriving at xjþ1=2 at time t with velocity u: ð1Þ

f0

      1 1 1 ð1Þ 1  H ut  xt2 ut  xt2 ; u  xt ¼ g l0 ð1 þ a1l ut  xt 2 2 2 2      1 1 ð1Þ H ut  xt 2 þ g r0 1 þ a1r ut  xt2 2 2 ð1Þ

ð1Þ

’ g l0 ð1  a1l utÞHðuÞ þ g r0 ð1  a1r utÞð1  HðuÞÞ; 2

where all the t terms are omitted. According to the relations between ð1Þ g l0 ð1Þ g r0

¼

ð1Þ g l ð1

l

 2k ðu 

U 1l Þ

¼

g ð1Þ r ð1

r

U 1r Þ

 2k ðu 

ð1Þ g l0 ;

ð19Þ

ð1Þ g r0

and

ð1Þ gl ;

ð1Þ gr

discussed above

xl tÞ; xr tÞ:

Substituting the above relation into Eq. (19),also omit the t2 terms, we have: ð1Þ

f0

  1 ð1Þ ut  xt2 ; u  xt ¼ g l ð1  2kl ðu  U 1l Þxl t  a1l utÞHðuÞ þ g rð1Þ ð1  2kr ðu  U 1r Þxr t  a1r utÞð1  HðuÞÞ: 2

ð20Þ

From the above equation, we can observe that the external force xl;r contributes to the variation of the particle distribution function on the first order of time t, which is on the the same order as the contributions from terms a1l and a1r . ð1Þ ð1Þ After the initial state f0 is constructed, the equilibrium state g 0 located at interface can be determined through the compatibility condition:

ZZ

þ1 1

ð1Þ Wð1Þ a g 0 du dn1 ¼

ZZ

ð1Þ Wð1Þ a f ðutÞdu dn1 ¼

ZZ u>0

ð1Þ Wð1Þ a g l0 du dn1 þ

ZZ u<0

ð1Þ Wð1Þ a g r0 du dn1 :

The equation above is equivalent to:

0

1

0

1

q10 q1l hu0 i>0 þ q1r hu0 i<0 C B 1 1C B q1l hu1 i>0 þ q1r hu1 i<0 A; @ q0 U 0 A ¼ @  1 2  2 2 1 1 2 1 1 q hu þ n i þ q hu þ n i q0 E0 >0 <0 r l 2

ð21Þ

where

q1l h   i>0 ¼ q1r h   i<0 ¼

Z u>0

Z

u<0

Z Z

ð1Þ

ð. . .Þg l0 du dn1 ; ð1Þ

ð. . .Þg r0 du dn1 :

The results of the moment of the Maxwellian distribution function can be found in [23]. According to Eq. (21), the conserð1Þ ð1Þ vative variables at the interface, which is denoted as W 0 , can be constructed. With the definition, the equilibrium states g 0 can be determined In order to determine the slope a1l and a1r , we make use of the relationship between distribution function and macroscopic variables:

Z Z

ð1Þ

1 ð1Þ Wð1Þ a al g 0 du dn1 ¼

1 ð1Þ Wð1Þ a ar g 0 du dn1 ¼

ð1Þ

2ðW 0  W j Þ h ð1Þ 2ðW jþ1

 h

The slope a1l can be computed from:

;

ð1Þ W0 Þ

:

7525

L. Pan et al. / Journal of Computational Physics 231 (2012) 7518–7536

0 11 a1l ð1Þ ð1Þ 1 2ðW 0  W j Þ B 1C ¼ M 10 ab @ a2l A; ð1Þ h q0 a13l where the matrix M 10 ab ¼

1

qð1Þ 0

R

ð22Þ

ð1Þ ð1Þ 1 Wð1Þ a Wb g 0 du dn1 . And ar can be computed similarly. ð1Þ

In the above equilibrium state g 0 , the external force effect has not been taken into consideration. As discussed above, we ð1Þ consider the relationship Eq. (15) between equilibrium state g ð1Þ ðx0 ; t0 ; u0 ; n1 Þ and g 0 ðx0 ; t 0 ; u0 ; n1 Þ. The integration of g ð1Þ ðx0 ; t 0 ; u0 ; n1 Þ in Eq. (16) becomes:

Z

t

0

0

g ð1Þ ðx0 ; t 0 ; u0 ; n1 Þeðtt Þ=s dt ¼

Z

0

t

0

ð1Þ

0

g 0 ½1  uðt  t 0 Þða1r ð1  HðuÞÞ þ a1l HðuÞÞ  2k0 ðu  U 10 Þx0 ðt  t0 Þeðtt Þ=s dt :

0

ð23Þ

Substituting Eqs. (20) and (23) into Eq. (16), we can obtain the distribution function of the solid phase at the cell interface: ð1Þ

ð1Þ

ð1Þ

f ð1Þ ðxjþ1=2 ; t; u; n1 Þ ¼ c0 g 0 þ c1 ½ða1r ð1  HðuÞÞ þ a1l HðuÞÞu  2k0 x0 ðu  U 0 Þg 0 þ c2 ðg l HðuÞ þ g rð1Þ ð1  HðuÞÞÞ ð1Þ

1 r r 1 þ c3 ðg l ð2kl ðu  U 1l Þxl þ a1l uÞHðuÞ þ g ð1Þ r ð2k ðu  U r Þx þ ar uÞð1  HðuÞÞÞ:

ð24Þ

where

c0 ¼ 1  et=s ; c1 ¼ sð1 þ et=s Þ þ tet=s ; c2 ¼ et=s ; c3 ¼ tet=s : After the distribution function for the solid phase f ð1Þ is determined, we need to construct the the distribution function for the gas phase f ð2Þ . Different from the solid phase, we need not to consider the effect of the external force, but the effect of the solid phase to the gas phase need to be taken into consideration. Eq. (7) can be rearranged as: ð2Þ

ft

ð2Þ

þ uf x ¼

sxfuð1Þ þ g ð2Þ  f ð2Þ : s

Denote

F1 ¼

1

s

Z

t

0

sf ð1Þ ðxjþ1=2 ; t0 ; u; n1 Þeðtt Þ=s dt0

ð25Þ

0

and

F2 ¼

1

s

Z 0

t

0

0

ð2Þ

g ð2Þ ðx0 ; t0 ; u; n2 Þeðtt Þ=s dt þ et=s f0 ðxjþ1=2  utÞ;

ð26Þ

where x0 ¼ xjþ1=2  uðt  t 0 Þ. The integral solution of Eq. (3.2) can be written as:

f ð2Þ ðxjþ1=2 ; t; u; nÞ ¼ F 2 ðx0 ; t 0 ; u; n2 Þ þ x0 ðF 1 Þu ðx0 ; t0 ; u; n1 Þ;

ð27Þ ð2Þ

where the x0 is the external force at the interface x ¼ xjþ1=2 . The initial distribution function f0 g ð2Þ around the interface can be constructed as: ð2Þ

and The equilibrium state

ð2Þ

2 f0 ¼ g l ð1 þ a2l ðx  xjþ1=2 ÞÞð1  Hðx  xjþ1=2 ÞÞ þ g ð2Þ r ð1 þ ar ðx  xjþ1=2 ÞÞHðx  xjþ1=2 Þ;

Fig. 1. These two pictures are taken from [2], which are the solid density and velocity of test 1.

ð28Þ

7526

L. Pan et al. / Journal of Computational Physics 231 (2012) 7518–7536 ð2Þ

g ð2Þ ¼ g 0 ð1 þ a2r ðx  xjþ1=2 Þð1  Hðx  xjþ1=2 ÞÞ þ a2l ðx  xjþ1=2 ÞHðx  xjþ1=2 ÞÞ;

ð29Þ

ð2Þ ð2Þ ð2Þ where g l ; g l and g 0 are local Maxwellian 2 2 2 2 interface, al ; ar , ar ; al are the corresponding

distribution functions located to the left, to the right and in the middle of a cell slopes. The dependence of a2l ; a2r ; a2r , a2l on the particle velocities is also obtained from the Taylor expansion of a Maxwellian and have the following form:

Fig. 2. These two pictures are taken from [8], which are the solid density and velocity of test 1.

0.315

second order first order exact

2.01

0.31

2

0.305

velocity

density

2.02

1.99

0.3

1.98

0.295

1.97

0.29

1.96

0.2

0.4

0.6

0.285

0.8

second order first order exact

0.2

0.4

x

0.6

0.8

x

0.8

12 second order first order exact

volume fraction

presure

10

8

second order first order exact

0.7

0.6

0.5

0.4 6

0.3 0.2

0.4

0.6

0.8

0.2

x

0.4

0.6

x Fig. 3. The numerical results of the first case for solid phase.

0.8

7527

L. Pan et al. / Journal of Computational Physics 231 (2012) 7518–7536

 1 2 u þ n22 ; 2  1 a2r ¼ a21r þ a22r u þ a23r u2 þ n22 ; 2  1 a2l ¼ a21l þ a22l u þ a23l u2 þ n22 ; 2  1 a2r ¼ a21r þ a22r u þ a23r u2 þ n22 ; 2 a2l ¼ a21l þ a22l u þ a23l

where all above coefficients are local constants. ð2Þ ð2Þ The equilibrium states g l and g r are obtained from macroscopic point-wise values at the cell interface:

Z

ð2Þ ð2Þ Wð2Þ a g l du dn2 ¼ W j ðxjþ1=2 Þ;

Z

ð2Þ ð2Þ Wð2Þ a g r du dn2 ¼ W jþ1 ðxjþ1=2 Þ:

ð30Þ

Then parameters ða21r ; a22r ; a23r ÞT and ða21r ; a22r ; a23r ÞT can be determined using the same method as the solid phase condition. The ð2Þ equilibrium state g 0 can be determined with the compatibility condition Eq. (8), we have:

ZZ

þ1 1

ð2Þ Wð2Þ a g 0 du dn2 ¼

ZZ

u>0

ð2Þ Wð2Þ a g l du dn2 þ

ZZ

u<0

ð2Þ Wð2Þ a g r du dn2 :

1

2.8

second order first order exact

2.6

velocity

density

0.8

0.6

0.4

second order first order exact

2.4

2.2

0.2 0.2

0.4

0.6

2

0.8

0.2

0.4

x

0.6

0.8

0.6

0.8

x

1 second order first order exact

0.8

presure

second order firsr order exact 0.6

0.4

volumefraction

0.8

0.6

0.4

0.2

0.2 0.2

0.4

0.6

0.8

0.2

x

0.4

x Fig. 4. The numerical results of the first case for gaseous phase.

7528

L. Pan et al. / Journal of Computational Physics 231 (2012) 7518–7536

The parameters ða21r ; a22r ; a23r ÞT and ða21r ; a22r ; a23r ÞT in g 0 can also be obtained by the same method discussed in the solid phase. All the parameter in f ð2Þ Eq. (28) have been determined, then F 2 can be written as ð2Þ

ð2Þ

ð2Þ

ð2Þ

2 2 F 2 ¼ c0 g 0 þ c1 ½ða2r ð1  HðuÞÞ þ a2l HðuÞÞug 0 þ c3 ðg ð2Þ r ð1  ar utÞð1  HðuÞÞ þ g l ð1  al utÞHðuÞÞ:

Substituting Eq. (24) into Eq. (25), we have: ð1Þ

ð1Þ

ð1Þ

ð1Þ

F 1 ¼ c00 g 0 þ c01 g 0 ½ða1r ð1  HðuÞÞ þ a1l HðuÞÞu  2k0 x0 ðu  U 10 Þ þ c02 ðg l HðuÞ þ g rð1Þ ð1  HðuÞÞÞ þ c03 ðg l ð2kl ðu  U 1l Þxl 1 r r 1 þ a1l uÞHðuÞ þ g ð1Þ r ð2k ðu  U r Þx þ ar uÞð1  HðuÞÞÞ;

where

c00 ¼ sð1  et=s Þ  tet=s ; 1 2

c01 ¼ s2 ð1 þ et=s Þ þ stet=s þ t2 et=s ; c02 ¼ tet=s ; 1 2

c03 ¼  t2 et=s :

1

1.4

1

velocity

density

1.2

second order first order exact

0.8

0.6

second order first order exact

0.8 0.6

0.4 0.4 0.2

0.2 0.2

0.4

0.6

0

0.8

0.2

0.4

x

0.6

0.8

x

1 0.8 second order first order exact

0.6

0.4

second order first order exact

volumefraction

presure

0.8

0.6

0.4

0.2

0.2

0.2

0.4

0.6

0.8

1

0.2

x

0.4

0.6

x Fig. 5. The numerical results of the second case for solid phase.

0.8

1

7529

L. Pan et al. / Journal of Computational Physics 231 (2012) 7518–7536

For simplicity, the corresponding slope terms of f ð1Þ and the t2 terms in Eq. (25) are omitted, we have: ð1Þ

ð1Þ

ð1Þ

F 1 ¼ c00 g 0 þ c001 g 0 ð2k0 x0 ðu  U 10 ÞÞ þ c02 ðg l HðuÞ þ g rð1Þ ð1  HðuÞÞÞ; where

c001 ¼ s2 ð1 þ et=s Þ þ stet=s : Substituting F 1 ; F 2 into Eq. (27), we obtain the gaseous distribution function at the cell interface. We have constructed the distribution functions of the solid and gaseous phase at the cell interfaces. The numerical fluxes can be computed by taking moments of the distribution functions: ðiÞ

F jþ1=2 ðtÞ ¼

Z

0 B u@

1

1 1 ðu2 2

C ðiÞ u Af ðxjþ1=2 ; t; u; nÞdu dni : 2 þ ni Þ

ð31Þ

ð2Þ

Especially, substituting fjþ1=2 ðtÞ into Eq. (31), we have: ð2Þ

F jþ1=2 ðtÞ ¼

Z

ð2Þ uWð2Þ a f du dn ¼

Z

0 uWð2Þ a ðF 2 þ x ðF 1 Þu Þdu dn2 ¼

Z

Z

0 uWð2Þ a F 2 du dn2  x

ðuWð1Þ a Þu F 1 du dn1 ;

R R ð1Þ where the integral ðuWð2Þ a Þu F 1 du dn2 was approximately considered as ðuWa Þu F 1 du dn1 .

1

1.4

1

velocity

density

1.2

second order first order exact

0.8

0.6

second order first order exact

0.8 0.6

0.4 0.4 0.2

0.2

0.2

0.4

0.6

0

0.8

0.2

0.4

x 1

0.8

1

second order first order exact

1

0.6

0.4

0.2

second order first order exact

0.8

volume fraction

0.8

presure

0.6

x

0.6

0.4

0.2

0.2

0.4

0.6

0.8

1

0.2

x

0.4

0.6

x Fig. 6. The numerical results of the second case for gaseous phase.

0.8

1

7530

L. Pan et al. / Journal of Computational Physics 231 (2012) 7518–7536

By integrating the above equations for a whole time step, we can obtain the numerical fluxes at the cell interfaces. 3.3. The non-conservative terms In this section, the source terms in Eqs. (13) and (14) can be treated explicitly. The update of flow variables is based on the following equations: nþ1;ð1Þ

Wj

n;ð1Þ

¼ Wj

n;ð2Þ

¼ Wj

þ þ

1 Dx 1 Dx

Z

t nþ1

Z

ð1Þ

ð1Þ

n;ð1Þ

; Wj

ð2Þ

ð2Þ

n;ð1Þ

;

ðF j1=2 ðtÞ  F jþ1=2 ðtÞÞdt þ DtSj

tn t nþ1

ðF j1=2 ðtÞ  F jþ1=2 ðtÞÞdt  DtSj

tn

ðiÞ

where F jþ1=2 ðtÞ is given by Eq. (31) and

@/ n;ð1Þ s

@x

j

n;ð1Þ

in Sj

nþ1;ð2Þ

has the following form:

 n;ð1Þ @/s ¼ Lððð/s Þj  ð/s Þj1 Þ=h; ðð/s Þjþ1  ð/s Þj Þ=hÞ; @x j where Lð; Þ is the Van Leer limiter defined in Eq. (3.2).

second order first order exact

2

second order first order exact

1

density

velocity

1.5

1

0.5

0

0.5

-0.5

0.2

0.4

0.6

0.8

0.2

0.4

0.6

6

0.2

5.5

second order first order exact

5

second order first order excat

0.18

4.5

volume fraction

4

presure

0.8

x

x

3.5 3 2.5 2 1.5

0.16

0.14

0.12

1 0.5 0.2

0.4

0.6

0.8

0.1

0.2

x

0.4

0.6

x Fig. 7. The numerical results of the third case the solid phase.

0.8

7531

L. Pan et al. / Journal of Computational Physics 231 (2012) 7518–7536

4. Numerical experiments In this section, we will present some numerical tests which have been tested in many papers [2,8,19,32]. These cases are presented to validate our algorithm for the numerical simulation of the two-phase flow calculations for the BN model. For all the test cases, the ideal gas EOS is used:

ea ¼

pa

qa ðca  1Þ

where ca ¼ 1:4; a ¼ s; g. For the ideal gas EOS, k ¼ q=2p. And for the stiffened gas EOS, if we denote e p a ¼ pa þ pa and consider the relationship ka ¼ qa =2 e p a , the expressions for the matrix M 1l ab and corresponding slopes are the same with the expressions of the ideal gas EOS. In the flux step, the collision time s at the interface is defined as:





!



pl  pr  pl  pr  g s s ¼ C 1 Dt þ C 2  gl Dt  þ  sl r r pg þ pg 

ps þ ps

where pai ; a ¼ s; g and i ¼ l; r represent the pressure of each phase at each side of the interface. In the test cases, C 1 = 0.05 and C 2 ¼ 25 are usually used. According to the numerical experiment, we observe that numerical results are not sensitive to the choices of the values of C 1 ; C 2 .

second order first order excat

1.1

second order first order exact

1.5

1 1

velocity

density

0.9

0.8

0.7

0.5

0

0.6 -0.5 0.5 0.2

0.4

0.6

0.8

0.2

0.4

x

0.6

0.8

x second order first order exact

3

0.95 second order first order exact

0.9

volumefraction

presure

2.5

2

0.8

1.5

1

0.85

0.2

0.4

0.6

0.8

0.75

0.2

x

0.4

0.6

x Fig. 8. The numerical results of the third case for gaseous phase.

0.8

7532

L. Pan et al. / Journal of Computational Physics 231 (2012) 7518–7536

The Van Leer limiter is used in the scheme for the reconstruction of conservative variables. CFL number equals to 0.15. The contours are the density, velocity, pressure, mass fraction of the solid and gaseous phase. In Test 1–Test 4, the meshes are n ¼ 400. In Test 5, we test the numerical scheme with n ¼ 400 and n ¼ 4000 meshes. Test 1: Single solid contact. Consider the following Riemann problem for the BN model: Phase

wL

rL

uL

pL

wR

qR

uR

pR

Solid Gas

0.8 0.2

2 1

0.3 2

5 1

0.3 0.7

2 0.1941

0.3 2.8011

12.8567 0.1

The solution of this problem consists of the single solid contact, propagating with the velocity 0.3. Numerical results of both solid and gaseous phase for first order scheme, second order scheme and the exact solution at t ¼ 0:1 are presented in Figs. 3 and 4 respectively. In Figs. 3 and 4, it shows that the numerical solution exhibits spurious oscillations near the solid contact, which are transported downstream in the course of time and could not be eliminated with improvement of the meshes. And it is caused by the presence of the nonconservative product. Compared with the first order scheme, we could observe that the numerical results was improved in the contour of the density and velocity of the solid phase. This problem was also tested in [2,8] by using the exact Riemann solvers and the numerical results are presented in Figs. 1 and 2. It could be noticed that the numerical solution of our second order scheme approximately equals to the exact solution except near the solid contact and the downstream in the contour of the density and velocity. From the relative error of average value in

0

2 BGK 1st order KFVS 1st order exact

1.8

1.6

BGK 1st order KFVS 1st order exact

velocity

density

-0.1

1.4

-0.2

1.2

1

-0.3 0.2

0.4

0.6

0.2

0.8

0.4

x

x

0.6

0.8

2 1.9

0.8

BGK 1st order KFVS 1st order exact

1.8

volume fraction

presure

BGK 1st order KFVS 1st order exact

0.7

1.7 1.6 1.5 1.4

0.6

0.5

1.3 1.2

0.4

1.1 1

0.2

0.4

0.6

0.8

0.2

0.4

0.6

x

x Fig. 9. The numerical results of the fourth case for the solid phase.

0.8

7533

L. Pan et al. / Journal of Computational Physics 231 (2012) 7518–7536

domains [0.35, 0.5] and [0.6, 0.75], it could be observed that our scheme performs well in this test. However, the maximum of overshot and undershot of our scheme at the solid contact is less accurate than the best-performing one in [8], which is associated with the treatment of non-conservative product. Test 2: We consider the following Riemann problem, and the initial data is taken as follows: Phase

wL

rL

uL

pL

wR

qR

uR

pR

Solid Gas

0.2 0.8

1 1

0.75 0.75

1 1

0.7 0.3

0.1 0.1

0 0

0.1 0.1

The problem was first tested in [21] and it is an extended version of a test in [22]. Numerical results of both solid and gaseous phase for first order scheme, second order scheme and the exact solution at t ¼ 0:15 are presented in Figs. 5 and 6 respectively. The solution, for each phase, contains a right shock wave, a contact discontinuity and a left sonic rarefaction wave. The correct resolution of the sonic point is very important for the numerical scheme. Compared with the method in [22], our current scheme also do not experience difficulty in resolving the sonic rarefaction wave. The method presented in [22] is an HLLC-type Riemann solver. Different from the HLLC-type Riemann solver for the Euler equations [22], they need to solve a nonlinear system which is based on the relation across the solid contact containing four equations. As a result, because of the nonlinear system, the iterative process is needed although it is a kind of approximate Riemann solver. In our scheme, numerical fluxes are obtained by taken moments of the distribution functions at the interface, so iterative process is avoided in our scheme, which is an advantage of our method. 0 BGK 1st order KFVS 1st order exact

1.4

-0.1

velocity

density

1.2

1

BGK 1st order KFVS 1st order exact

-0.2

0.8 -0.3 0.6 0.2

0.4

0.6

0.8

0.2

0.4

x 0.7

2

BGK 1st order KFVS 1st order exact

0.8

BGK 1st order KFVS 1st order exact

0.6

volumefraction

presure

1.8

1.6

1.4

1.2

1

0.6

x

0.5

0.4

0.3

0.2

0.2

0.4

0.6

0.8

0.1

0.2

0.4

x

0.6

x Fig. 10. The numerical results of the fourth case for gaseous phase.

0.8

7534

L. Pan et al. / Journal of Computational Physics 231 (2012) 7518–7536

Test 3: Coinciding shocks and rarefaction waves. We consider following problem: Phase

wL

rL

uL

pL

wR

qR

uR

pR

Solid Gas

0.1 0.9

0.2068 0.5806

1.4166 1.5833

0.0416 1.375

0.2 0.8

2.2263 0.4890

0.9366 0.7011

6 0.986

A particular issue on the Riemann problem for the BN model is that several waves can coincide with each other. However, the gas parameters do not change across the solid waves. Thus, the numerical solution across these waves should be independent of the presence of the wave of the other phase. The solution of this problem consists of two coinciding 1-shock waves for the gaseous and solid phase, and the gaseous 3-shock inside of the solid 3-rarefaction. Numerical results of of density, velocity, pressure and fraction of both solid and gaseous phase for first order scheme, second order scheme and the exact solution at t ¼ 0:15 are presented in Figs. 7 and 8 respectively. It could be observed that the gaseous phase presents a good agreements with the exact solution. The undershoots in solid density and velocity are observed, and most previous schemes exhibits the same behavior [2]. Test 4: The comparison of the BGK scheme and the KFVS scheme. We consider the following Riemann problem: Phase

wL

rL

uL

pL

wR

qR

uR

pR

Solid Gas

0.4 0.6

1 0.5

0 0

1 1

0.8 0.2

2 1.5

0 0

2 2

2

-0.2

n=4000 n=400

n=4000 n=400

1.5

uelo

den

-0.4

-0.6

-0.8 1 -1

0.2

0.4

0.6

0.2

0.8

0.4

0.6

0.8

x

x 0.6

n=4000 n=400

n=4000 n=400

volumefraction

pre

4

3

0.4

0.2

2

0.2

0.4

0.6

0.8

0

0.2

x

0.4

0.6

x Fig. 11. The numerical results of the fifth case for solid phase.

0.8

7535

L. Pan et al. / Journal of Computational Physics 231 (2012) 7518–7536

In the test, we make the comparison between the first order of the current BGK scheme and KFVS scheme [32]. The first order numerical results for BGK scheme, KFVS scheme and the exact solution at t ¼ 0:1 are presented in Figs. 9 and 10 respectively. It could be observed that, except the oscillation in the solid velocity, BGK scheme could approximate the exact solutions better. Especially, in the contact discontinuity in Fig. 10, the smeared area shrinks. Both of BGK scheme and KFVS scheme are based on the gas kinetic theory. As analyzed in [23], in the gas evolution stage, KFVS scheme is based on the collisionless model and BGK scheme is based on the generalized BGK model which contains the collision terms. In KFVS scheme, the artificial collision term with the artificial collision time, which equals to the time step, helps the KFVS scheme to capture these waves. Test 5: Gas shock approaches solid contact. We consider the following problem: Phase

wL

rL

uL

pL

wR

qR

uR

pR

Solid Gas

0.5 0.5

2.1917 6.3311

usL 0.789

3 1

0.1 0.9

0.6333 0.4141

1.1421 0.6741

2.5011 0.0291

The ill-posedness of the Riemann problem for the Baer–Nunziato model was discussed in [2,8]. There are some initial conditions for which the Riemann problem admits more than an unique entropy-satisfying solutions or for which it does not admit a solution at all. It is interesting to mention the above initial condition where usL ¼ 0:995, it correspond to those of Test 4 in[2], for which the authors observed differences between numerical results and exact solution. Meanwhile, accord-

6

-0.2 n=4000 n=400

5

-0.4

uelo

4

den

n=4000 n=400

-0.3

3

-0.5

-0.6

2

-0.7

1 0.2

0.4

0.6

0.8

0.2

0.4

0.8

1

1

0.9

n=4000 n=400

volume fraction

0.8

presure

0.6

x

x

0.6

0.4

0.2

0.8

n=4000 n=400

0.7

0.6

0.5

0.2

0.4

0.6

0.8

0.4

0.2

x

0.4

0.6

x Fig. 12. The numerical results of the fifth case for gaseous phase.

0.8

7536

L. Pan et al. / Journal of Computational Physics 231 (2012) 7518–7536

ing to the discussion in [8], when 2:46 < usL < 1:51, there is no solution to the Riemann problem. As a matter of fact, usL ¼ 0:995 is a case of existence of two solutions and the two solution of this problem were presented in [8] . In Figs. 11 and 12, we present the second order results at t ¼ 0:15 for the BGK scheme with the n ¼ 4000 and n ¼ 400 meshes. With the improvement of the mesh, it could be observed that the numerical results of the current scheme could converge to a steady solution which corresponds to one of the solutions of the Riemann problem. The ill-posedness of the Riemann problem for Baer–Nunziato model is an interesting problem which needs to be research further in the future. 5. Conclusion In this paper, we have developed a new gas kinetic scheme for the homogeneous Baer–Nunziato two-phase flow model containing non-conservative terms. In the view of gas kinetic theory, a generalized BGK model which matches with the Baer– Nunziato model is proposed. Based on this model, we construct the the numerical scheme. Numerical tests are presented in this paper and validate the current approach for the simulation of the BN model. Compared with the exact and approximate Riemann solver, in our scheme, numerical fluxes can be obtained by taking moments of the distribution functions, and nonconservative terms are explicitly introduced into the construction of numerical fluxes. In consequence, the iterative process of exact solutions is not needed in this method. Compared with the collisionless KFVS scheme, the collisions of particles are taken into consideration. The ill-posedness of the Riemann problem for the Baer–Nunziato model was also discussed. In our scheme, the numerical results could converge to one of the solutions of this Riemann problem. It could be noticed that, in the generalized BGK model, we only consider the collisions between the same species so that we may lose some useful physical information. In the future work, we will present numerical method based on a new BGK model in which it takes the collisions both with particles of their species and the other species into consideration. References [1] R. Abgrall, How to prevent pressure oscillations in multicomponent flow calculations: a quasi conservative approach, J. Comput. Phys. 125 (1996) 150– 160. [2] N. Andrianov, G. Warnecke, The Riemann problem for the Baer–Nunziato two-phase flow model, J. Comput. Phys. 195 (2004) 434–464. [3] M.R. Baer, J.W. Nunziato, A two-phase mixture theory for the deflagration-to-detonation transition in reactive granular materials, Int. J. Multiphase Flow. 12 (1986) 861–889. [4] J.B. Bdzil, R. Menikoff, S.F. Son, A.K. Kapila, D.S. Stewart, Two-phase modelling of deflagration-to-detonation transition in granular materials: a critical examination of modelling issues, Phys. Fluids 11 (1999) 378–402. [5] P.L. Bhatnagar, E.P. Gross, M. Krook, A model for collision processes in gases I: small amplitude processes in charged and neutral one-component systems, Phys. Rev. 94 (1954) 511–525. [6] C. Cercignani, The Boltzmann Equation and its Applications, Springer-Verlag, 1988. [7] S. Chapman, T.G. Cowling, The Mathematical Theory of Non-Uniform Gases, third ed., Cambridge University Press, 1990. [8] V. Deledicque, M.V. Papalexandris, An exact Riemann solver for compressible two-phase flow models containing non-conservative products, J. Comput. Phys. 222 (2007) 217–245. [9] P. Embid, M. Baer, Mathematical analysis of a two-phase continuum mixture theory, Contin. Mech. Thermodyn. 4 (1992) 279–312. [10] S. Gavrilyuk, R. Saurel, A compressible multiphase model with microinertia, J. Comput. Phys. 175 (2002) 326–360. [11] S. Jiang, G.X. Ni, A c-model BGK scheme for compressible multifluids, Int. J. Numer. Methods Fluids 46 (2004) 163–182. [12] P.G. Lefloch, A.E. Tzavaras, Representation of weak limits and definition of nonconservative products, SIAM J. Math. Anal. 30 (1999) 1309–1342. [13] R. Saurel, O. Lemetayer, A multiphase model for compressible flows with interfaces, shocks, detonation waves and cavitation, J. Fluid Mech. 431 (2001) 239–271. [14] Q. Li, K. Xu, S. Fu, A high-order gas-kinetic Navier–Stokes flow solver, J. Comput. Phys. 229 (2010) 6715–6731. [15] Y.S. Lian, K. Xu, A gas-kinetic scheme for multimaterial flows and its application in chemical reactions, J. Comput. Phys. 163 (2000) 349–375. [16] G. Dal Maso, P.G. Lefloch, F. Murat, Definition and weak stability of nonconservative products, J. Math. Pures Appl. 74 (1995) 483–548. [17] M.V. Papalexandris, A two-phase model for compressible granular flows based on the theory of irreversible processes, J. Fluid Mech. 517 (2004) 103– 112. [18] R. Saurel, R. Abgrall, A multiphase Godunov method for compressible multifluid and multiphase flows, J. Comput. Phys. 150 (1999) 425–467. [19] D.W. Schwendeman, C.W. Wahle, A.K. Kapila, The Riemann problem and a high-resolution Godunov method for a model of compressible two-phase flow, J. Comput. Phys. 212 (2006) 490–526. [20] H.Z. Tang, K. Xu, A high-order gas-kinetic method for multidimentional ideal magnetohydrodynamics, J. Comput. Phys. 165 (2000) 69–88. [21] S.A. Tokareva, E.F. Toro, HLLC-Type Riemann solver for the Baer–Nunziato compressible two-phase flow, J. Comput. Phys. 229 (2010) 3573–3604. [22] E. Toro, Riemann Solvers and Numerical Methods for Fluid Dynamics, Springer, 1997. [23] K. Xu, Gas Kinetic Schemes for Unsteady Compressible Flow Simulations, Lecure Note Ser., vol. 1998-03, Von Karman Institute for Fluid Dynamics Lecture, 1998. [24] K. Xu, A gas-kinetic BGK scheme for the Navier–Stokes equations and its connection with artificial dissipation and Godunov method, J. Comput. Phys. 171 (2001) 289–335. [25] K. Xu, Numerical Hydrodynamics from Gas-Kinetic Theory, PhD. thesis, Columbia University, 1993. [26] K. Xu, BGK-based scheme for multicomponent flow calculations, J. Comput. Phys. 134 (1997) 122–133. [27] K. Xu, Discontinuous Galerkin BGK method for viscous flow equations: one-dimensional systems, SIAM J. Sci. Comput. 23 (2004) 1941–1963. [28] K. Xu, Gas-kinetic thoery based flux splitting method for ideal magnetohydrodynamics, J. Comput. Phys. 153 (1999) 344–375. [29] K. Xu, E. Josyula, Continuum formulation for non-equilibrium shock structure caculation, Commun. Comput. Phys. 1 (2006) 425–450. [30] K. Xu, K. Prendergast, Numerical Navier–Stokes solutions from gas-kinetic theory, J. Comput. Phys. 114 (1994) 9–17. [31] K. Xu, A well-balanced gas-kinetic scheme for the shallow-water equations with source terms, J. Comput. Phys. 178 (2002) 533–562. [32] B. She, G. Zhao, A Gas-kinetic scheme for compressible two-phase flow model containing non-conservative products, Chinese J. Comput. Phys. 29 (2012) 51–57.