DOI 10.1007/s00193-006-0037-y O R I G I NA L A RT I C L E
A volume-fraction based algorithm for hybrid barotropic and non-barotropic two-fluid flow problems
Keh-Ming Shyue
Received: 9 December 2004 / Accepted: 25 January 2006 / Published online: 13 July 2006
© Springer-Verlag 2006
Abstract The aim of this paper is to describe a simple Eulerian interface-capturing approach for the efficient numerical resolution of a hybrid barotropic and non-barotropic two-fluid flow problem in more than one space dimension. We use the compressible Euler equa- tions as a model system with the thermodynamic prop- erty of each of the barotropic and non-barotropic fluid components characterized by the Tait and Noble–Abel equations of state, respectively. The algorithm is based on a volume fraction formulation of the equations to- gether with an extended equation of state that is de- vised to give an approximate treatment for the mixture of more than one fluid component within a grid cell. A standard high-resolution wave propagation method is employed to solve the proposed two-fluid model with the dimensional-splitting technique incorporated in the method for multidimensional problems. Several numer- ical results are presented in one and two space dimen- sions that show the feasibility of the algorithm as applied to a reasonable class of practical problems without the occurrence of any spurious oscillation in the pressure near the smeared material interfaces. This includes, in particular, solutions for a study on the variation of the jet velocity with the incident shock pressure arising from the collapse of an air cavity in water under a shock wave.
Keywords Volume-fraction based algorithm· Tait equation of state· Noble–Abel equation of state · Gas-cavity collapse· Jetting
Communicated by K. Takayama.
K.-M. Shyue (
B
)Department of Mathematics, National Taiwan University, Taipei 106, Taiwan
e-mail: [email protected]
PACS 02.70.Bf· 46.15.-x · 47.11.+j · 47.55.Kf
1 Introduction
This paper is concerned with a simplified compressible two-phase flow problem where there is a material inter- face separating regions of two different barotropic and non-barotropic fluid components within an Nd≥ 1 spa- tial domain. In this problem, the flow regime of interest is assumed to be homogeneous with no jumps in the pres- sure and velocity (the normal component of it) across the interface, and the physical effects such as the viscos- ity, surface tension, and heat conduction are assumed to be small, and hence can be ignored in the problem formulation. With that, we use an Eulerian viewpoint of the model system, that on the barotropic part of the domain, the fluid is governed by the isentropic version of the compressible Euler equations as
∂
∂t
ρ ρui
+
Nd
j=1
∂
∂xj
ρuj
ρuiuj+ p(ρ)δij
= 0, (1)
for i= 1, 2, . . . , Nd, while on the non-barotropic part of the domain, the fluid is governed by the full set of the compressible Euler equations as
∂
∂t
ρ ρui
ρE
+Nd
j=1
∂
∂xj
ρuj
ρuiuj+ pδij
ρEuj+ puj
= 0. (2)
Here ρ, uj, p, E, and δij denote the density, the parti- cle velocity in the xj-direction, the pressure, the specific total energy, and the Kronecker delta, respectively.
To complete the model, in this work, the constitu- tive law for each of the barotropic and non-barotropic
components is taken to satisfy the Tait equation of state for compressible liquids (or called the Murnaghan equa- tion of state in the context of an elastic solid [25]), p(ρ) = (p0+ B)
ρ ρ0
γ
− B, (3)
and the Noble–Abel equation of state for real gases (or called the constant covolume equation of state [40]), p(ρ, e) =
γ − 1 1− bρ
ρe, (4)
in a respective manner. Note that in the above expres- sions e represents the internal energy, ρ0 and p0 the density and pressure at a reference state, respectively,B a pressure-like constant, b a finite size of volume occu- pied by molecules (0≤ b < 1/ρ), and γ the adiabatic constant (γ > 1), see [20,23,27,45] for a typical set of data of practical importance. As usual, we have E= e + Nd
j=1u2j/2. Sample applications of this two-fluid model are, for instance, to the simulation of a shock-induced collapse of gas cavities in liquid [7, 8] or to the simulation of an underwater explosion bubble [21, 44].
For an efficient numerical resolution of this hybrid barotropic and non-barotropic two-fluid flow problem, we want to use a generalization of the classical shock- capturing method designed originally for single compo- nent flows in an Eulerian framework. It is known in the literature that the principal difficulty in the usual mod- ification is the occurrence of spurious pressure oscilla- tions when two or more fluid components are present in a grid cell. In the case of a barotropic or non-barotropic multicomponent problem, there has been quite a few numerical methods developed for that matter, see [1–
3, 9, 14, 15, 26, 31, 43] and the references cited therein for more exposition. In the current interest of the problems consisting of separate barotropic and non-barotropic re- gions, there is, however, a relatively few attempts de- voted to the subject, see [5, 21, 37] for an example using a somewhat complicated Lagrangian-type approach as compared to the current Eulerian one.
The basic approach we take here is an extension of the previous work advocated by the author [33, 34, 36] in that a mixture equation of state is introduced first as a basis to the modeling of the mixing between two different barotropic and non-barotropic fluid components. Then with the help of a volume-fraction function, we define an extended equation of state so that the pressure of the fluid can be determined explicitly no matter what fluid component (pure or not) is within a grid cell (see Sect. 2.1). Having that, as before, we are able to derive a volume-fraction based on the model system that consists of the full Euler equations for the basic conserved vari- ables and an additional set of equations for the volume
fraction as well as the partial density of any one of the two fluid components in the problem. In our method, the latter equations are included in the algorithm primarily for an easy determination of the approximate location of the interface and also the computation of the pressure from the equation of state. It is important to note that these equations are put in a form so as to ensure a consis- tent modeling of the energy equation near the smeared interfaces and also the fulfillment of the mass equation in the other single component regions (see Sect. 2.2).
To find approximate solutions of our model system, we use a high-order Godunov method based on a wave- propagation formulation with the dimensional-splitting technique incorporated in the method for multidimen- sional problems. Numerical results to be presented in Sect. 4 show that this is a viable approach to a reason- able class of practical problems in one and two space dimensions, without producing any wrong oscillations in the pressure near the interfaces.
It should be emphasized that the methodology we have proposed here is by no means limited to the two- fluid flows with the above chosen equations of state.
Extension of the algorithm to problems involving more than two fluid components and more complicated equa- tions of state can be made in a straightforward manner by following the idea described in this paper and the five-equation model of Allaire et al. [3] and Massoni et al. [24]. Without going into the details for that, our goal is to establish the basic solution strategy and val- idate its use via some sample numerical experimenta- tions.
The format of this paper is as follows. In Sect. 2, we describe the basic mathematical equations in a fluid- mixture form that will be used later in a numerical method for constructing approximate solutions of our model hybrid barotropic and non-barotropic two-fluid problems. In Sect. 3, we give a brief review of the HLLC approximate Riemann solver and the high-resolution wave-propagation methods for solving problems in one and more than one space dimensions. Numerical results of some sample examples are shown in Sect. 4.
2 Mathematical formulation
2.1 Equations of state
To start out, we want to build a mixture equation of state that is prerequisite in our approach for modeling the numerical mixing between the barotropic and non- barotropic fluid components within a grid cell. The first step, we take here, is to assume that thermodynamically the mixture of the fluids characterized by (3) and (4)
behaves like a typical non-barotropic fluid. Then using the first and second laws of thermodynamics, we may rewrite (3) and (4) in terms of the specific entropy S and the specific volume V= 1/ρ as
p(V, S) = A(S) (p0+ B)
V0
V γ
− B, (5)
and
p(V, S) = A(S)p0
V0− b V− b
γ
, (6)
for each in turn, where A(S) = R exp [(S − S0)/CV] (cf. [10]). Here CV represents the specific heat at con- stant volume,R the universal gas constant, and S0the specific entropy at a reference state. As in the previ- ous work for multicomponent problems with a van der Waals equation of state [33], we assume further that all the fluid components under concern is in an adiabatic equilibrium with the same entropy, and accordingly we propose to define an equation of state for the fluid mix- ture that combines (5) with (6) to be of the form p(V, S) = A(S) (p0+ B)
V0− b V− b
γ
− B. (7)
Notice that, with the use of this modified Noble–Abel equation of state (7), we recover not only the case (6) whenB → 0, but also the case (5) and so (3) as we set S= S0in (5), when b→ 0. In the non-limiting case, how- ever, where none of the material-dependent quantities B and b is close to zero, (7) does give a way to the repre- sentation of the cases in between, that is to the mixing of the present barotropic and non-barotropic components.
It should be mentioned that as a practical matter for many fluid flow problems away from cavitation (which is the situation studied here), rather than using (7) for the fundamental thermodynamic description of the fluid mixture, it is customary to formulate the equation in a way with the often-used variables in compressible hydrodynamics:ρ and e instead. Again, by employing the basic thermodynamic principles, this can be done quite easily, leading to
p(ρ, e) =
γ − 1 1− bρ
(ρe − B) − B. (8)
As to the computation of the fluid temperature T which is important in a wide variety of applications, we may use one of the formulas,
p(ρ, T) = ρRT 1− bρ − B, e(ρ, T) = RT
γ − 1+B ρ, for realization.
Now, for an easy identification of the type of fluid component within a grid cell, a volume-fraction func- tionα is introduced here for that purpose. For example, when grid cells contain only the barotropic component we may setα = 1, and so when grid cells contain only the non-barotropic component we setα = 0. In case there are some cells cut by the interfaces where α ∈ (0, 1), we then have both of these components occupied by the volume fractionsα and 1 − α for each separately. With this definition ofα, the pressure of our two-fluid flow problem in all the fluid-component scenarios within a grid cell can be determined straightforwardly by
p=
(p0+ B)
ρρ0
γ
− B if α = 1
γ −1 1−bρ
(ρe − B) − B if α = 1, (9)
provided all the variables appeared there are defined and known a priori. Finally, it should be remarked that in this work, the thermodynamical description of the materials of interest is limited by the stability require- ment that the speed of sound of the fluid belongs to a set of real numbers.
2.2 Equations of motion
With the proposed constitutive law (9), it is clear that in regions whereα = 0 or α = 1, there is no problem in using (1) and (2), respectively, for the complete descrip- tion of the behavior of the underlying single-compo- nent flow. In regions whereα ∈ (0, 1), as the fluid mix- ture is modeled in a non-barotropic manner, it is quite common to apply (2) as a model system that describes the motion of the mixtures of the basic conserved vari- ables. Besides that, as in our previous work on numeri- cal methods for compressible multicomponent problems (cf. [32, 33, 35, 36]), we also bring in a set of effective equations for the problem-dependent material quanti- ties so that the pressure can be computed easily from the equation of state.
2.2.1γ -based effective equations
To derive the aforementioned effective equations for the mixture of material quantities in the current hybrid barotropic and non-barotropic flow application, we start with an interface-only problem as usual where both the pressure and each component of the particle velocities are constants in the domain, while the other variables such as the density and the material quantities are having jumps across some interfaces. Then, from (2), it is easy
to obtain equations for the time-dependent behavior of the density and total internal energy as
∂ρ
∂t +
Nd
j=1
uj∂ρ
∂xj = 0,
∂
∂t(ρe) +
Nd
j=1
uj ∂
∂xj (ρe) = 0,
in a respective manner. Now, by inserting the modified Noble–Abel equation of state (8) into the latter one, we find an alternative description of the energy equation
∂
∂t
1− bρ
γ − 1p+γ − bρ γ − 1 B
+
Nd
j=1
uj ∂
∂xj
1− bρ
γ − 1p+γ − bρ γ − 1 B
= 0 (10)
that is in relation to not only the pressure, but also the density and material quantities:γ , b, and B.
In our algorithm, to maintain the pressure in equilib- rium as it should be for this interface-only problem, we split (10) into the following two equations for the fluid mixtures(1 − bρ)/(γ − 1) and (γ − bρ)B/(γ − 1) as
∂
∂t
1− bρ γ − 1
+
Nd
j=1
uj ∂
∂xj
1− bρ γ − 1
= 0, (11)
∂
∂t
γ − bρ γ − 1 B
+
Nd
j=1
uj ∂
∂xj
γ − bρ γ − 1 B
= 0, (12)
respectively. We emphasize that in order to have the correct pressure equilibrium in (10), these are the two key equations that should be satisfied and approximated consistently (when the problem is solved numerically, see Sect. 3). On the other hand, as a practical matter, it is obvious that in addition to (11) and (12), we need to impose one more additional condition so as to have enough equations for the determination of all the three material quantities in (8). In our approach, see [33], this is done by simply breaking (11) into the following two parts:
∂
∂t
1
γ − 1
+
Nd
j=1
uj ∂
∂xj
1
γ − 1
= 0, (13)
∂
∂t
bρ γ − 1
+
Nd
j=1
uj ∂
∂xj
bρ γ − 1
= 0. (14)
Having done so, we arrive at a system of three equations (12)–(14) for the variables(γ − bρ)B/(γ − 1), 1/(γ − 1), and bρ/(γ − 1), in the order mentioned. Combining them to (2) yields a model system that is fundamental in
our algorithm for describing the behavior of the numer- ical mixing between the barotropic and non-barotropic components near the interface. With that, there is no difficulty to compute the pressure based on the equa- tion of state,
p=
ρE − Nd
j=1(ρuj)2
2ρ −γ − bρ γ − 1 B
1− bρ
γ − 1
.
Up to this point, our discussion stresses only on an approach that is capable of keeping the pressure in equilibrium for a model interface-only problem. Since in practice we are interested in shock wave problems as well, we should take the equations, i.e., (12), (13), and (14), in a form such that γ , b, and B remain un- changed across both shocks and rarefaction waves. In this regard, it is known that with γ governed by (13), there is no problem to do so (cf. [1, 32]). For b andB, however, due to the appearance of the density term in (12) and (14), it turns out that in a time when such a situation occurs, for consistency with the mass conser- vation law of the fluid mixture they should be modified by
∂
∂t
γ − bρ γ − 1 B
+
Nd
j=1
∂
∂xj
γ − bρ γ − 1 Buj
= γ B γ − 1
Nd
j=1
∂uj
∂xj
, (15)
and
∂
∂t
bρ γ − 1
+
Nd
j=1
∂
∂xj
bρ γ − 1uj
= 0, (16)
respectively, so that the mass-conserving property of the solution in the single component region can be acquired also (cf. [33]). We note that for convenience, we call the set of equations: (13), (15), and (16), aγ -based effective equations for the mixture of the material quantities of the modified Noble–Abel equation of state to be distinct from the other one presented below.
2.2.2α-based effective equations
Before proceeding further, we comment that to find the initial fluid mixtures, 1/(γ − 1), (γ − bρ)B/(γ − 1), and bρ/(γ − 1), that is necessary when we initialize the data for multicomponent flow computations, we use the equation of state (8), where written as a function of the volume fractionαι, forι = 1, 2, and α1+ α2 = 1, it reads
1− bρ
γ − 1p+γ − bρ
γ − 1 B = ρe =
2 ι=1
αιριeι
=
2 ι=1
αι
1− bιρι
γι− 1 pι+γι− bιρι γι− 1 Bι
.
Here the subscript “ι” denotes the state variable of fluid componentι. By taking a similar approach as employed in Sect. 2.2.1 for the derivation of theγ -based effective equation it comes out readily as a splitting of the above expression into the relations:
1 γ − 1 =
2 ι=1
αι
γι− 1, γ − bρ γ − 1 B =
2 ι=1
αιγι− bιρι γι− 1 Bι, bρ
γ − 1 =
2 ι=1
αι bιρι
γι− 1, (17)
where in the process of splitting the terms the pressure p is chosen to fulfill the condition
1− bρ γ − 1
p=
2 ι=1
αι
1− bιρι γι− 1
pι. (18)
Notice that when each of the partial pressures is in an equilibrium state within a grid cell, in conjunction with the first and third parts of (17), the pressure p obtained from (18) would remain in the same equilibrium as well, i.e., p= pι, forι = 1, 2.
Now if the above volume-fraction notion of the states 1/(γ − 1), (γ − bρ)B/(γ − 1), and bρ/(γ − 1) are being employed in theγ -based effective equations, we are able to rewrite them straightforwardly into a componentwise form as
∂
∂t
αι γι− 1
+
Nd
j=1
uj ∂
∂xj
αι γι− 1
= 0, (19)
∂
∂t
αιγι− bιρι γι− 1 Bι
+
Nd
j=1
∂
∂xj
αιγι− bιρι γι− 1 Bιuj
= αι γιBι
γι− 1
Nd
j=1
∂uj
∂xj
, (20)
∂
∂t
αι bιρι
γι− 1
+
Nd
j=1
∂
∂xj
αι bιρι
γι− 1uj
= 0, (21)
forι = 1, 2. Then based on the fact that all the material quantitiesγι, bι, andBιwill be kept as a constant in each phase of the domain at all time, from (19), it is easy to find the transport equation for the volume fractionαιas
∂αι
∂t +
Nd
j=1
uj∂αι
∂xj = 0, (22)
whereas from (21), we find the mass conservation law for the fluid componentι as
∂
∂t(ριαι) +
Nd
j=1
∂
∂xj
ριαιuj
= 0. (23)
Furthermore, in the case of (20), after some simple alge- braic manipulations, it can be expressed by the form
γιBι γι− 1
∂αι
∂t +
Nd
j=1
uj∂αι
∂xj
− bιBι γι− 1
∂∂t(ριαι) +
Nd
j=1
∂
∂xj
ριαιuj
= 0,
yielding naturally the same set of governing equations as before, i.e., (22) and (23). It is apparent that if the solutions ofαι andριαι are known from the equations, we may therefore compute 1/(γ −1), (γ −bρ)B/(γ −1), and bρ/(γ − 1) directly according to (17). Thus, instead of using the γ -based effective equations, it is a viable alternate to use the α-based equations: (22) and (23), for the motion of the mixture of the material quantities of the problem.
2.2.3 Complete model system
It is important to mention that to make up a complete model system that is capable of dealing with all the fluid component cases, α = 1 (barotropic phase), α = 0 (non-barotropic phase), or α ∈ (0, 1) (barotropic and non-barotropic coexistence phase), we have to know the approximate location of the interfaces so that the correct equations of motion as well as the equation of state can be employed to each part of the domain, from the current time to the next. Here, asα1+ α2 = 1, it is clear that if we chooseα1 = α and so α2 = 1 − α, the two transport equations in (22) for each ofα1andα2can be combined, without affecting anything, to a single one forα as
∂α
∂t +
Nd
j=1
uj∂α
∂xj = 0, (24)
leading to the evolution equation we use in the algo- rithm for that matter. Note that, in devising a fluid-mix- ture type algorithm for multicomponent problems, one common practice is to consider the mixture of the total density,ρ = ρ1α1+ ρ2α2, as one of the basic variables in the proposed model system (cf. [32, 33, 35]). When this
is the case, it should be more sensible to include only one of the equations in (23) together with the last two relations in (17) than the use of both (15) and (16) for the variables(γ − bρ)B/(γ − 1) and bρ/(γ − 1).
Putting all things together, with the equation of state (9), the model equations with which we propose to solve the present two-fluid flow problems in more than one space dimension take the form
∂ρ
∂t +
Nd
j=1
∂
∂xj
ρuj
= 0
∂
∂t(ρui) +
Nd
j=1
∂
∂xj(ρuiuj+ pδij) = 0
∂
∂t(ρE) +
Nd
j=1
∂
∂xj
ρEuj+ puj
= 0 if α = 1, (25)
∂
∂t(ρ1α) +
Nd
j=1
∂
∂xj
ρ1αuj
= 0
∂α
∂t +
Nd
j=1
uj∂α
∂xj = 0
for i= 1, 2, . . . , Nd. Note that whenα = 1 we have a sys- tem of Nd+ 4 equations for the motion of both the pure non-barotropic phase and the fluid mixture. Clearly, the first Nd+ 2 of them are simply the basic conservation of the mass, momenta (Nd of them), and total energy, while the last two equations are (23) and (24) that are used to find the variables 1/(γ − 1), (γ − bρ)B/(γ − 1), and bρ/(γ − 1) from the relations in (17), and also the approximate location of the interface (the volume frac- tion functionα does it). On the other instance, when α = 1, we just have the first Nd+ 1 equations governing the single-component barotropic flow as usual. With a system expressing in this way, there is no problem to compute all the state variables of interest, including the pressure from the equation of state. The initialization of the state variables in (25) for fluid-mixture cells can be made in a standard way as described in [32, 33] for numerical simulations.
It should be remarked that for an easy extension of the proposed model to more complex equations of state we have taken (25) in a form that is analogous to the five- equation model advocated by Allaire et al. [3] or Mas- soni et al. [24] for compressible multicomponent flow problems with general non-barotropic equation of state.
For any given Nd, if the state variables of the flow are all in the region of the thermodynamic stability (this is the case we are interested in here), it is not difficult to show that (25) is a hyperbolic system in the sense that any
linear combination of the matrices Aj, j= 1, 2, . . . , Nd, appearing in the quasi-linear form of the equations
∂q
∂t +
Nd
j=1
Aj(q)∂q
∂xj = 0 (26)
has real eigenvalues and a complete set of eigenvectors.
As an example, we consider the three-dimensional case Nd = 3, and then have the state vector q in (26) defined by
q= (ρ, ρu1, ρu2, ρu3, ρE, ρ1α, α)T, and the matrices Aj, for j= 1, 2, 3, defined by
A1=
0 1 0 0 0 0 0
K− u21 2u1+ η2 η3 η4 ϕ χ
−u1u2 u2 u1 0 0 0 0
−u1u3 u3 0 u1 0 0 0
u1(K − H) H + u1η2u1η3u1η4u1( + 1) u1ϕ u1χ
−u1ρ1α/ρ ρ1α/ρ 0 0 0 u1 0
0 0 0 0 0 0 u1
,
A2=
0 0 1 0 0 0 0
−u1u2 u2 u1 0 0 0 0
K− u22 η2 2u2+ η3 η4 ϕ χ
−u2u3 0 u3 u2 0 0 0
u2(K − H) u2η2H+ u2η3u2η4u2( + 1) u2ϕ u2χ
−u2ρ1α/ρ 0 ρ1α/ρ 0 0 u2 0
0 0 0 0 0 0 u2
,
A3=
0 0 0 1 0 0 0
−u1u3 u3 0 u1 0 0 0
−u2u3 0 u3 u2 0 0 0
K− u23 η2 η3 2u3+ η4 ϕ χ u3(K − H) u3η2u3η3H+ u3η4u3( + 1) u3ϕ u3χ
−u3ρ1α/ρ 0 0 ρ1α/ρ 0 u3 0
0 0 0 0 0 0 u3
.
With that, the eigenvalues and the corresponding eigen- vectors of the matrices are:
for matrix A1:
A1 = diag
λ(1)1 , λ(1)2 , . . . , λ(1)7
= diag(u1− c, u1, u1+ c, u1, . . . , u1), RA1 =
r(1)1 , r(1)2 ,. . . , r(1)7
=
1 1 1 0 0 0 0
u1− c u1 u1+ c 0 0 0 0
u2 u2 u2 1 0 0 0
u3 u3 u3 0 1 0 0
H− u1c K/ H + u1c u2u3−ϕ/ −χ/ ρ1α/ρ 0 ρ1α/ρ 0 0 1 0
0 0 0 0 0 0 1
,
for matrix A2: A2 = diag
λ(2)1 , λ(2)2 , . . . , λ(2)7
= diag(u2− c, u2, u2+ c, u2, . . . , u2), RA2 =
r(2)1 , r(2)2 ,. . . , r(2)7
=
1 1 1 0 0 0 0
u1 u1 u1 1 0 0 0
u2− c u2 u2+ c 0 0 0 0
u3 u3 u3 0 1 0 0
H− u2c K/ H + u2c u1u3−ϕ/ −χ/ ρ1α/ρ 0 ρ1α/ρ 0 0 1 0
0 0 0 0 0 0 1
,
and
for matrix A3: A3 = diag
λ(3)1 , λ(3)2 , . . . , λ(3)7
= diag(u3− c, u3, u3+ c, u3, . . . , u3), RA3 =
r(3)1 , r(3)2 ,. . . , r(3)7
=
1 1 1 0 0 0 0
u1 u1 u1 1 0 0 0
u2 u2 u2 0 1 0 0
u3− c u3 u3+ c 0 0 0 0 H− u3c K/ H + u3c u1u2−ϕ/ −χ/
ρ1α/ρ 0 ρ1α/ρ 0 0 1 0
0 0 0 0 0 0 1
;
Ajr(j)k = λ(j)k r(j)k , j = 1, 2, 3, and k = 1, 2, . . . , 7. Here c =
[γ /(1 − bρ)][(p + B)/ρ] is the speed of sound of the fluid, and the other notations appeared in the above formulae are set by = (γ − 1)/(1 − bρ), K = (3
j=1u2j/2) + ζ , H = E + p/ρ, ηi = −ui, for i = 1, 2, 3,ϕ = [bB/(γ − 1)] − p[b/(γ − 1)]/ , and χ = −[γ B/(γ − 1)] + p[1/(γ − 1)]/ , where ζ =
−b2/(γ2− 1)[(B2) + (p/ )] and [κ] = κ1− κ2. It should be noted that since the volume fraction func- tionα is governed by a linear transport equation of the form (24) and its solution can only have jumps across the material interfaces, but not across other waves such as shock and rarefaction, we easily find the usual expres- sion of the Rankine–Hugoniot jump conditions across shock waves for this two-fluid model (25) (cf. [33]).
Finally, for the ease of the latter reference, it is cus- tomary to write (25) into a more compact expression by
∂q
∂t +
Nd
j=1
fj
∂
∂xj
, q
= 0, (27)
where fjis taken as the vector-value function of the form fj=
∂
∂xj
ρuj
, ∂
∂xj
ρu1uj+ pδ1j
, . . . , ∂
∂xj
ρuNduj
+pδNdj
, ∂
∂xj
ρEuj+ puj
, ∂
∂xj
ρ1αuj
, uj∂α
∂xj
T
, for j= 1, 2, . . . , Nd.
3 Numerical approximations
To find approximate solutions of our two-fluid model presented in Sect. 2, we use a standard high-resolution wave propagation method developed by LeVeque [17, 18] for general hyperbolic systems of partial differen- tial equations. This method is a different form of the fluctuation-and-signal scheme of Roe [29, 30] in that we solve one-dimensional Riemann problems at each cell interface, and use the resulting waves (i.e., discontinu- ities moving at constant speeds) to update the solutions in neighboring grid cells. To achieve second-order accu- rate on smooth solutions, and sharp and monotone pro- files on discontinuous solutions, we introduce slopes and limiters to the method as in many other high-resolution schemes for conservation laws [11, 19, 39].
3.1 HLL-type approximate Riemann solvers
As a preliminary, we begin by reviewing the construction of the approximate Riemann problem solutions in one space dimension which is one of the major elements in our numerical algorithm. If we consider the case Nd= 3 as an example, to simplify the notation, rather than using xjfor the spatial variable and ujfor the particle veloc- ity, for j = 1, 2, 3, we take the often-employed symbols (x, y, z) and (u, v, w) for that matter instead. Then the problem to be solved is to find the solution of (25) in the direction normal to one of the xy, yz, and xz planes, with piecewise constant data qLand qRto the left and right of the cell interface. It should be noted that in the applications concerned here, we have chosen the data qL and qRwell enough so that the solution of the Riemann problem would consist of genuinely non-linear waves such as shock and rarefaction, and linearly degenerate wave such as contact discontinuity; there is no vacuum region occurring in the solution. Without loss of gen- erality, in the following we only look at the Riemann problem in the direction normal to the yz plane.
There are a wide variety of approaches proposed in the literature for determining an approximate solution of the Riemann problem, see [41], for example. Here we are interested in a simple variant of the Riemann solver
based on the work of Harten et al. [13] for hyperbolic systems of conservation laws,
∂q
∂t + ∂
∂xf(q) = 0, (28)
where q∈ lRnis the vector of conserved variables for a system of n equations and f the flux function. Recall that in the original version of the HLL solver, the solution of the Riemann problem is assumed to be composed of two discontinuities propagating at constant speeds λL and λR to the left and right,λL = λR, separating three constant states in the space-time domain. If we assume further thatλL and λR are known a priori by some simple estimates based on the local information of the wave speeds, then by using the integral form of the conservation laws over a sufficiently large control vol- ume[−M, M] × [0, T], for some positive M and T ∈ lR, it is easy to find the constant state in the middle region, denoted by qm, as the average of the exact solution over the interval[λLT,λRT] at time T,
qm = 1
(λR− λL)T
λRT
λLT
q(x, T) dx
=λRqR− λLqL− f (qR) + f (qL) λR− λL
, (29)
where f(qι) is the flux evaluated at the state qι, forι = L or R.
In spite of the fact that the above 2-wave HLL solver has been applied quite successfully to many numerical methods for solving problems governed by (28), it is known, however, that the numerical result obtained by using this solver is too diffusive for contact discontinu- ities. In addition to that, because of the lack of infor- mation on the structure of the material interfaces, it is not feasible at all to the more general multicomponent problems like the one considered in this article. Nev- ertheless, as recommended by Toro et al. [6, 42], this 2-wave Riemann solver can be improved quite easily by introducing an additional middle wave of speed um in the solution structure for modeling the speed of contact discontinuity, yielding a 3-wave HLL (or called HLLC) solver. Note that if we have had the first two compo- nents of qm computed from (29), we may simply set um = q(2)m /q(1)m (cf. [31, 36]), where q(i)m is the ith com- ponent of the vector qm, see [41] for the various other possible choices.
With that, our goal next is to find the constant states qmL and qmR in the regions mL and mR to the left and right of the middle wave, respectively. We do this by making use of the integral form of the conserva- tion laws again, but now applied over a control volume [−M, umT − ] × [0, T] for the state qmL, and over a
control volume[umT+ , M] × [0, T] for the state qmR, 0 < 1. When the aforementioned procedure is taken effect to the current two-fluid model (25) in the prescribed direction, it is not difficult to show that the result is
q(1)mι = fm(1)ι λι− um
, q(2)mι = umq(1)mι, q(3)mι = vιq(1)mι,
q(4)mι = wιq(1)mι, q(5)mι = fmι(5)+ um
umfmι(1)− fmι(2) λι− um
, (30)
q(6)mι = fm(6)ι λι− um
, q(7)mι = q(7)ι ,
where fmι(i) = λιq(i)ι − f(i)(qι) represents the ith compo- nent of the vector fmι, forι = L or R. It is interesting to remark that q(i)mL and q(i)mR satisfy the basic consistency condition of the integral form of the conservation laws,
um− λL
λR− λL
q(i)Lm+
λR− um
λR− λL
q(i)mR = q(i)m,
for i= 1, 2, . . . , 6. By taking λL= min (uL− cL, uR− cR) andλR= max (uL+ cL, uL+ cL), for instance, we may then set the speed of the three moving discontinuities by λ1 = λL,λ2 = um,λ3 = λR, and the jumps across each of them by
W1= qmL− qL, W2= qmR− qmL, W3= qR− qmR, using the result (30). As usual, wave propagation meth- ods (to be described below) are based on using these propagating discontinuity to update the cell averages in the cells neighboring each interface.
3.2 Wave propagation methods 3.2.1 One-dimensional case
To review the basic idea of our underlying integration scheme, it is instructive to look at the simplest Nd = 1 case of our model system,
∂q
∂t + f1
∂
∂x, q
= 0,
with q and f1defined as in (25) and (27). For simplicity, we describe the method on a uniform grid with fixed mesh spacingx in the x-direction, and use a standard finite-volume formulation in which the discretized solu- tion Qnj approximates the cell average over the grid cell [xj, xj+1] at time tn,
Qnj ≈ 1
x
xj+1
xj
q(x, tn) dx.
The time step from the current time tnto the next tn+1 is denoted byt.
In this numerical setup, a first order accurate version of the method in wave-propagation form is a Godunov- type scheme that can be written as
Qnj+1= Qnj − t
x
mw
m=1
λ−mWm
n j+1+
λ+mWm
n
j . (31)
Hereλm and Wm are solutions of the mth wave fam- ily, for m = 1, 2, . . . , mw, obtained from solving the Riemann problems at cell interfaces xj and xj+1 with a properly chosen solver, and λ− = min (λ, 0), λ+ = max(λ, 0). It is known that method (31) belongs to a class of upwind schemes and is stable when the typical CFL (Courant–Friedrichs–Lewy) condition is satisfied (cf. [11, 19]). Moreover, it is not difficult to show that the method is quasi-conservative in the sense that when applying the method to (25) not only the conservation laws but also the transport equation is approximated in a consistent manner by the method, see [33] for the details.
To improve the accuracy of (31) to high resolution, it is a standard approach by first introducing correction waves in a piecewise-linear form with zero mean value, and then propagating each wave over the time stept to update the cell averages it overlaps. Without providing the details in this place (cf. [19] for example), with the corrections, (31) is modified by
Qn+1j := Qn+1j −1 2
mw
m=1
|µm| (1 − |µm|) Wm
n
j+1
−
|µm| (1 − |µm|) Wm
n
j , (32)
whereµm= λmt/x, and Wmis a limited value ofWm
obtained by comparingWmwith the correspondingWm
from the neighboring Riemann problem to the left (if λm> 0) or to the right (if λm< 0).
Note that when utilizing the wave-form representa- tion of the Riemann problem solution, it is common practice to perform the limiting procedure over each component of the wave via a limiter function (eg, by using the minmod function(θ) = max(0, min(1, θ)) or some others as discussed in [39]), and set
Wmj(i)= θmj(i)
Wmj(i) with θmj(i)=WmJ(i) Wmj(i), J=
j− 1 if λmj≥ 0
j+ 1 if λmj< 0, (33)
where Wmj(i) is the ith component of Wmj. While the above approach works in a satisfactory manner for many single component problems, it was shown in [36] that,
for general multicomponent problems, to ensure a consistent approximation of the total energy near the interfaces the third limited component of the 2-wave, W2j(3), should be modified so that the desired pressure equilibrium can be maintained there.
Having this in mind, following the previous work done in [36], we make a new definition of W2j(3)as
W2j(3) :=
W¯2J
W¯2j
W¯2j−
bB γ − 1
W2j(4)
+
γ B γ − 1
W2j(5),
where ¯W2ι = W2ι(3)+ [bB/(γ − 1)]W2ι(4)− [γ B/(γ − 1)]W2ι(5), forι = j and J (ie, j − 1 or j + 1 depending on the propagating direction ofλ2j);[κ] = κ1− κ2. With this revision of the limited 2-wave on the total energy, it is not difficult to show that for a model interface- only problem (see Sect. 2.2) we again have the required pressure equilibrium without introducing any spurious oscillations in the approximate solution. Besides that we would typically obtain a better resolution of the result as compared to the first order result, see [36] for details.
3.2.2 Multidimensional case
To extend the one-dimensional wave propagation method to more space dimensions, here we take a simple dimensional-splitting approach in which a multidimen- sional problem is split into a sequence of one-dimen- sional problems. Consider the three-dimensional case Nd = 3, for example. The hybrid barotropic and non- barotropic two-fluid flow problem modeled by (27), i.e.,
∂q
∂t + f1
∂
∂x, q
+ f2
∂
∂y, q
+ f3
∂
∂z, q
= 0,
can be split into
x-sweeps: ∂q
∂t + f1
∂
∂x, q
= 0, (34a)
y-sweeps: ∂q
∂t + f2
∂
∂y, q
= 0, (34b)
z-sweeps: ∂q
∂t + f3
∂
∂z, q
= 0. (34c)
Assume a uniform Cartesian grid with fixed mesh spac- ingx, y, and z in the x-, y-, and z-direction, respec- tively. In a finite-volume formulation of the solution, the value Qnijk would approximate the cell average of the solution over the(i, j, k)th grid cell at time tn,
