A fluid-mixture type algorithm for barotropic two-fluid flow problems
Keh-Ming Shyue
Department of Mathematics, National Taiwan University, Taipei 106 Taiwan, ROC Received 28 January 2004; received in revised form 28 April 2004; accepted 5 May 2004
Available online 20 July 2004
Abstract
Our goal is to present a simple interface-capturing approach for barotropic two-fluid flow problems in more than one space dimension. We use the compressible Euler equations in isentropic form as a model system with the thermodynamic property of each fluid component characterized by the Tait equation of state. The algorithm uses a non-isentropic form of the Tait equation of state as a basis to the modeling of the numerically induced mixing between two different barotropic fluid components within a grid cell. Similar to our previous work for multicomponent problems, see [J.
Comput. Phys. 171 (2001) 678] and references cited therein, we introduce a mixture type of the model system that consists of the full Euler equations for the basic conserved variables and an additional set of evolution equations for the problem-dependent material quantities and also the approximate location of the interfaces. A standard high-resolution method based on a wave-propagation formulation is employed to solve the proposed model system with the dimen- sional-splitting technique incorporated in the method for multidimensional problems. Several numerical results are presented in one, two, and three space dimensions that show the feasibility of the method as applied to a reasonable class of practical problems without introducing any spurious oscillations in the pressure near the smeared material interfaces.
Ó 2004 Elsevier Inc. All rights reserved.
AMS: 65M06; 76L05; 76M20; 76T05
Keywords: Barotropic flow; Tait equation of state; Stiffened gas equation of state; Fluid-mixture model; Wave propagation method
1. Introduction
In this paper, we are concerned with a model barotropic two-fluid flow problem where the flow regime of interest is assumed to be homogeneous with no jumps in the pressure and velocity (the normal component of it) across a material interface in an NdP1 spatial domain. In the problem formulation, we use the is- entropic version of the compressible Euler equations of the form
E-mail addresses:shyue@math.ntu.edu.tw,shyue@jacobi.math.ntu.edu.tw (K.-M. Shyue).
0021-9991/$ - see front matter Ó 2004 Elsevier Inc. All rights reserved.
doi:10.1016/j.jcp.2004.05.003
www.elsevier.com/locate/jcp
o ot
q qui
þXNd
j¼1
o oxj
quj
quiujþ pðqÞdij
¼ 0 for i¼ 1; 2; . . . Nd ð1Þ
for the basic equations of motion of each fluid component, where q, uj, p, and dij denote the density, the particle velocity in the xj-direction, the pressure, and the Kronecker delta function, respectively (cf. [9]). To complete the model, the constitutive law for the thermodynamic property of the fluid is taken to satisfy the Tait equation of state for compressible liquids (or called the Murnaghan equation of state in the context of an elastic solid [34]),
pðqÞ ¼ pð 0þ BÞ q q0
c
B: ð2Þ
Here p0and q0are the reference pressure and density, respectively, B is a pressure-like constant, and c is a dimensionless parameter, see [31] for a typical set of material-dependent quantities of practical importance.
Representative applications of this two-fluid model are, for instance, to the simulation of the propagation of shock waves in a water–human tissue media in Extracorporeal Shock Wave Lithotripsy [35], or in a water–silicone oil media in the semiconductor industry [57].
To solve this barotropic flow problem numerically, we want to use a generalization of the classical shock-capturing method designed originally for single component flows. For non-barotropic multicom- ponent problems, it is known in the literature that the principal difficulty in the usual extension is the occurrence of spurious pressure oscillations when two or more fluid components are present in a uniform Cartesian grid cell, see the sample numerical methods proposed in [1–3,12,15,18,21,30,43] for more ex- position. For the current application of barotropic flows, there is, however, a relatively few methods de- veloped for the matter, see the recent work of Koren et al. [23], and van Brummelen and Koren [54] for a concise survey.
The approach we take here is to first introduce a non-isentropic form of the Tait equation of state as a basis to the modeling of the mixing between two different barotropic fluid components. Then with the help of a volume-fraction function, we define a hybrid 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 Section 2).
Having that, as in the previous work [45–48], we are able to derive a mixture type of the model system that consists of the full Euler equations for the basic conserved variables and an additional set of evolution equations for the problem-dependent material quantities and also the approximate location of the inter- faces. In our approach, the latter equations are included in the algorithm primarily for an easy computation of the pressure from the equation of state, and are put in a form so as to ensure a consistent modeling of the energy equation near the smeared interfaces, and also the fulfillment of the mass equation in the other single component regions. A standard high-resolution method based on a wave-propagation formulation is employed to solve the proposed model system with the dimensional-splitting technique incorporated in the method for multidimensional problems. Numerical results to be presented in Section 6 show that this is a viable approach to a reasonable class of practical problems without producing any wrong oscillations in the pressure near the interfaces.
It should be emphasized that the approach we have proposed here is by no means limited to the two-fluid flows with a Tait equation of state. Extension of the algorithm to other barotropic flow problems with the more general pressure law as occurred in cases with some solid materials (cf. [32,33]), for example, and also to problems involving more than two-flow components, can be made in a similar manner by following the idea described in this paper. Without going into the details for that here, we want to focus our attention to the establishment of the basic solution strategy and validate its use via numerical experimentation of some sample problems.
This paper is organized as follows. In Section 2, we propose a new equation of state for the modeling of the mixing between two different barotropic flows within a grid cell. In Section 3, we describe the model equations that is easy to use for numerical approximation of barotropic two-fluid problems. The con- struction of the approximate solution for the one-dimensional Riemann problem of the model is discussed in Section 4, and a brief review of numerical methods based on wave propagation is given in Section 5. In Section 6, we present a sample test of results for problems in one, two, and three space dimensions.
2. Equations of state
One of the key elements in our fluid-mixture approach for multicomponent problems is the introduction of a hybrid version of the equation of state that is necessary in the algorithm for modeling the numerical mixing between two different fluid components within a grid cell. In the current instance of barotropic flows, the basic idea of the proposed method is simple. We view the mixture of two different barotropic fluids, where each of them is described by the same Tait equation of state (2) but possibly with a different set of material-dependent quantities: c, B, and q0, as a non-barotropic fluid, and use an entropy-dependent equation of state for describing the thermodynamic behavior of the fluid mixture instead. By using the first and second laws of thermodynamics (cf. [13,59]) together with (2), it is easy to derive such an equation of state for the pressure p in terms of the mixture density q and the specific entropy S as:
pðq; SÞ ¼ AðSÞ pð 0þ BÞ q q0
c
B; ð3Þ
where AðSÞ ¼ exp½ðS S0Þ=CV and CVdenotes the specific heat at constant volume. Clearly, (3) reduces to the barotropic flow case (2), when the change of the entropy DS¼ S S0! 0. In the non-limiting case, however, when DS is not close to zero, (3) does give a way to the representation of the cases in between, that is to the mixing of two different barotropic fluids. Numerical experiences indicate that when DS is not large, i.e., entropy variation of the fluids is somewhat small, (3) is a good model to use for a reasonable class of practical problems, see Section 6.
It is worthwhile to mention that when we rewrite (3) in terms of the often-used variables in gas dynamics, qand the internal energy e, it takes the form
pðq; eÞ ¼ ðc 1Þq e
þB
q0
cB ð4Þ
and when we rewrite (3) in terms of q and the temperature T , it becomes pðq; T Þ ¼ qRT B;
where R is the universal gas constant. Note that in the fluid dynamical literature an equation of state of the form (4) is typically called the stiffened gas or Tammann equation of state (cf. [14,17]). It is also easy to show that the internal energy of the purely barotropic fluid can be computed by (4) but with the use of (2) directly for the pressure term.
For convenience in the later development, we use a volume-fraction function Y to indicate the type of fluid component within a grid cell (this is a standard way to do in many two-phase flow solver [42,56,58]).
For instance, when grid cells contain only fluid-component 1 we may set Y ¼ 1, and so when grid cells contain only fluid-component 2 we set Y ¼ 0. In case there are some cells cut by the interfaces where Y 2 ð0; 1Þ, we then have both the fluid-components 1 and 2 occupied by the volume fractions Y and 1 Y , respectively. With this definition of Y , the pressure of a barotropic two-fluid flow problem in all the fluid- component scenarios within a grid cell can be determined straightforwardly by:
p¼
p01þ B1
ð Þ qq
01
c1
B1 if Y ¼ 1 ðbarotropic fluid 1Þ;
p02þ B2
ð Þ qq
02
c2
B2 if Y ¼ 0 ðbarotropic fluid 2Þ;
ðc 1Þq e þBq
0
cB if Y 2 ð0; 1Þ ðmixture of two barotropic fluidsÞ 8>
>>
<
>>
>:
ð5Þ
provided that 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 requirement that the speed of sound of the fluid belongs to a set of real numbers. Certainly, it is both interesting and important to include the cavitation and phase-transition effects to the present constitutive model in a region where the pressure drops to a critical value, but this subject matter is a very difficult one, and is beyond the scope of this paper, see [4,5,11,44,55] and references cited therein for some possible models and approaches proposed in the literature.
3. Equations of motion
We now discuss the model equations that will be used in our numerical methods for constructing ap- proximate solutions of barotropic two-fluid problems. To begin, clearly, in regions where Y ¼ 0 or Y ¼ 1, it is enough to use (1) and (2) for the complete description of the behavior of the underlying single-component barotropic flow. In regions where Y 2 ð0; 1Þ, however, since the mixture of two barotropic fluid components is modeled in a non-isentropic manner, see Section 2, the original isentropic Eq. (1) are not sufficient to this instance, and so some supplementary equations need to be considered to provide further information of the fluid mixture.
It is apparent that, because the thermodynamic behavior of the fluid mixture depends on the entropy as we have postulated, the conservation law for the total energy of the form
o
otðqEÞ þXNd
j¼1
o oxj
qEuj
þ puj
¼ 0 ð6Þ
should be incorporated in the model system, where E¼ e þPNd
j¼1u2j=2 is the specific total energy. In ad- dition to that, as in our previous work on numerical methods for compressible multicomponent problems (cf. [45,46,48]), we also introduce a set of effective equations for the problem-dependent material quantities so that the pressure can be computed easily from the equation of state.
3.1. c-based effective equations
To derive the aforementioned effective equations for the mixture of material quantities in the present stiffened gas case, it is usual to start with an interface-only problem where both the pressure and each component of the particle velocities are constant in the domain, while the other variables such as the density and the material quantities are having jumps across some interfaces. In this instance, from (1) and (6), it is easy to obtain equations for the time-dependent behavior of the density and total internal energy as
oq ot þXNd
j¼1
ujoq oxj
¼ 0;
o
otðqeÞ þXNd
j¼1
uj o oxj
ðqeÞ ¼ 0;
in a respective manner. By inserting the equation of state (4) into the latter one, we find an alternative description of the energy equation
o ot
pþ cB c 1
B q0q
þXNd
j¼1
uj o oxj
pþ cB c 1
B q0q
¼ 0 ð7Þ
that is in relation to not only the pressure, but also the material quantities: c, B, and q0.
In our algorithm, to maintain the pressure in equilibrium as it should be for our model interface-only problem, we split (7) into the following two equations for the fluid mixtures of 1=ðc 1Þ and cB=ðc 1Þ Bq=q0as
o ot
1 c 1
þXNd
j¼1
uj o oxj
1 c 1
¼ 0; ð8aÞ
o ot
cB c 1
B q0q
þXNd
j¼1
uj o oxj
cB c 1
B q0q
¼ 0; ð8bÞ
respectively. We emphasize that in order to have the correct pressure equilibrium in (7), these are the two key equations that should be satisfied and approximated consistently (when the problem is solved numerically, see Section 5). On the other hand, as a practical matter, it is obvious that, in addition to (8a) and (8b), we need to impose one more additional condition so as to have enough equations for the three material quantities: c, B, and q0. In our approach (cf. [46,48]), this is done by simply splitting (8b) into the following two parts:
o ot
cB c 1
þXNd
j¼1
uj o oxj
cB c 1
¼ 0; ð8cÞ
o ot
B q0q
þXNd
j¼1
uj o oxj
B q0q
¼ 0: ð8dÞ
Having done so, we arrive at a system of three equations. (8a), (8c), and (8d) for the variables 1=ðc 1Þ, cB=ðc 1Þ, and Bq=q0, respectively. Combining them to (1) and (6) yields a model system that is fun- damental in our algorithm for describing the behavior of the numerical mixing between two barotropic fluids near the interface. With that, there is no difficulty to compute the pressure from (4) as
p¼ qE
"
PNd
j¼1ðqujÞ2
2q þB
q0q cB c 1
# 1
c 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., (8a), (8c), and (8d), in a form such that c, B, and q0remain unchanged across both shocks and rarefaction waves. In this regard, it is easy to see that with c and B governed by (8a) and (8c), respectively, there is no problem to do so (cf. [1,45]). For q0 or B=q0, however, due to the ap- pearance of the density term in (8d), it turns out that, in a time when such a situation occurs, for consistent with the mass conservation law of the fluid mixture the primitive form of (8d) should be modified by
o ot
B q0q
þXNd
j¼1
o oxj
B q0quj
¼ 0; ð8eÞ
so that the mass-conserving property of the solution in the single component region can be acquired also (cf. [46]). We note that, for convenience, we call the set of equations (8a), (8c), and (8e), a c-based effective
equations for the mixture of the material quantities of the stiffened gas equation of state to be distinct from the other one presented below.
3.2. Y -based effective equations
Before proceeding further, we note that to find the initial fluid mixtures, 1=ðc 1Þ, cB=ðc 1Þ, and Bq=q0, that is necessary when we initialize the data for multicomponent flow computations, we use the equation of state (4), where written as a function of the volume fraction Yi, for i¼ 1; 2, and Y1þ Y2¼ 1, it reads
pþ cB c 1 B
q0q¼ qe ¼X2
i¼1
Yiqiei¼X2
i¼1
Yi piþ ciBi
ci 1
Bi
q0iqi
: ð9aÞ
Here the subscript ‘‘i’’ denotes the state variable of fluid component i. By taking a similar approach as employed in Section 3.1 for the derivation of the c-based effective equation (cf. [45,46] also), it comes out easily a splitting of (9a) into the following set of relations:
1
c 1¼X2
i¼1
Yi
ci 1; cB c 1¼X2
i¼1
Yi
ciBi
ci 1; and B
q0q¼X2
i¼1
Yi
Bi
q0iqi; ð9bÞ
where in the process of splitting the terms the pressure p is chosen to satisfy the relation as p
c 1¼X2
i¼1
Yi pi
ci 1: ð9cÞ
With the first part of (9b), it is easy to see that when each of the partial pressures is in equilibrium within a grid cell, the pressure obtained from (9c) would remain in equilibrium also, i.e., p¼ pi, for i¼ 1; 2, see [37]
for a different way to derive the same result.
Now with the volume-fraction notion of the states 1=ðc 1Þ, cB=ðc 1Þ, and Bq=q0being defined by (9b), the c-based effective equations may be rewritten into the form:
o ot
X2
i¼1
Yi ci 1
! þXNd
j¼1
uj
o oxj
X2
i¼1
Yi ci 1
!
¼ 0; ð10aÞ
o ot
X2
i¼1
Yi ciBi
ci 1
! þXNd
j¼1
uj o oxj
X2
i¼1
Yi ciBi
ci 1
!
¼ 0; ð10bÞ
o ot
X2
i¼1
YiBi
q0iqi
! þXNd
j¼1
o oxj
X2
i¼1
YiBi
q0iqiuj
!
¼ 0: ð10cÞ
After some simple algebraic manipulations (cf. [45]), from both (10a) and (10b), we find the transport equation for each volume fraction Yi as
oYi
ot þXNd
j¼1
ujoYi
oxj
¼ 0; ð10dÞ
while from (10c) we find the mass conservation equation for each fluid component i as o
otðYiqiÞ þXNd
j¼1
o oxj
Yiqiuj
¼ 0; ð10eÞ
for i¼ 1; 2. Clearly, when the set of Yiand Yiqiare known from (10d) and (10e), we may therefore compute 1=ðc 1Þ, cB=ðc 1Þ, and Bq=q0directly from (9b). Thus, instead of using the c-based effective equations, it is a viable approach to use the Y -based equations (10d) and (10e), for the motion of the mixture of the material quantities of the problem.
3.3. Complete model system
Note that to constitute a complete model system that is capable of dealing with all the fluid phase cases, Y ¼ 0, Y ¼ 1, or Y 2 ð0; 1Þ, 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, since Y1þ Y2¼ 1, it is clear that if we choose Y1¼ Y and so Y2¼ 1 Y , the two transport equations in (10d) for each of Y1 and Y2 can be combined, without affecting anything, to a single one for Y as
oY ot þXNd
j¼1
uj
oY oxj
¼ 0; ð10fÞ
leading to the evolution equation we use in practice for that matter. It is important to mention that, in devising a fluid-mixture type algorithm for multicomponent problems, one common practice (cf. [45,46,48]) is to consider the mixture of the total density, q¼ Y1q1þ Y2q2, as one of the basic variables, but is not to use the mixtures of the separate fluid densities, Y1q1and Y2q2, as the basic variables, in the proposed model system. When this is the case, it should be more sensible to use (8e) as the governing equation for the time- evolution of the variable Bq=q0than the use of the last relation in (9b) together with one of the equations in (10e) for the determination of Bq=q0.
Putting all the things together, with the hybrid equation of state (5), the model equations with which we propose to solve the barotropic two-fluid flow problems in more than one space dimension take the form:
oq otþXNd
j¼1
o oxj
ðqujÞ ¼ 0;
o
otðquiÞ þXNd
j¼1
o oxj
ðquiujþ pdijÞ ¼ 0 for i¼ 1; 2; . . . ; Nd; o
otðqEÞ þXNd
j¼1
o oxj
qEuj
þ puj
¼ 0 if Y 2 ð0; 1Þ;
o ot
B q0q
þXNd
j¼1
o oxj
B q0quj
¼ 0 if Y 2 ð0; 1Þ;
oY ot þXNd
j¼1
ujoY oxj
¼ 0:
ð11Þ
Notice that when Y 2 ð0; 1Þ we have a system of Ndþ 4 equations for the motion of the mixing between two barotropic phases. Clearly, the first Ndþ 2 of them are simply the basic conservation of the mass, momenta (Ndof them), and total energy, while the last two equations are (8e) and (10f) that are introduced to model the fluid-mixing of Bq=q0 and the volume-fraction function Y (a quantity that is used in the model to identify the approximate location of the interface and also to find c and B according to the first two re- lations in (9b)). On the other instance, when Y ¼ 0 or Y ¼ 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 (11) for fluid-mixture cells can be made in a standard way as described in [45,46] for numerical simulations.
Note that, as before, the proposed model system (11) is not written in the full conservation form, but is rather a quasi-conservative system of equations. However, 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 (11) 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
oq otþXNd
j¼1
AjðqÞoq
oxj¼ 0 ð12aÞ
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 (12a) defined by
q¼ q; qu1;qu2;qu3;qE;B q0q; Y
T
; ð12bÞ
and the matrices Aj, for j¼ 1; 2; 3, defined by
A1¼
0 1 0 0 0 0 0
K u21 u1ð2 CÞ u2C u3C C C v
u1u2 u2 u1 0 0 0 0
u1u3 u3 0 u1 0 0 0
u1ðK H Þ H u21C u1u2C u1u3C u1ðC þ 1Þ u1C u1v
u1B=q0 B=q0 0 0 0 u1 0
0 0 0 0 0 0 u1
2 66 66 66 66 66 64
3 77 77 77 77 77 75
;
A2¼
0 0 1 0 0 0 0
u1u2 u2 u1 0 0 0 0
K u22 u1C u2ð2 CÞ u3C C C v
u2u3 0 u3 u2 0 0 0
u2ðK H Þ u1u2C H u22C u2u3C u2ðC þ 1Þ u2C u2v
u2B=q0 0 B=q0 0 0 u2 0
0 0 0 0 0 0 u2
2 66 66 66 66 66 64
3 77 77 77 77 77 75
;
A3¼
0 0 0 1 0 0 0
u1u3 u3 0 u1 0 0 0
u2u3 0 u3 u2 0 0 0
K u23 u1C u2C u3ð2 CÞ C C v u3ðK H Þ u1u3C u2u3C H u23C u3ðC þ 1Þ u3C u3v
u3B=q0 0 0 B=q0 0 u3 0
0 0 0 0 0 0 u3
2 66 66 66 66 66 64
3 77 77 77 77 77 75 :
With that, the eigenvalues and the corresponding eigenvectors of the matrices are:for matrix A1, KA1¼ diag k ð1Þ1 ;kð1Þ2 ; . . . ;kð1Þ7
¼ diagðu1 c; u1; u1þ c; u1; . . . ; u1Þ;
RA1¼ r1ð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=C Hþ u1c u2 u3 1 v=C
B=q0 0 B=q0 0 0 1 0
0 0 0 0 0 0 1
2 66 66 66 66 66 64
3 77 77 77 77 77 75
;
for matrix A2,
KA2¼ diag k ð2Þ1 ;kð2Þ2 ; . . . ;kð2Þ7
¼ diagðu2 c; u2; u2þ c; u2; . . . ; u2Þ;
RA1¼ r1ð2Þ; 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=C Hþ u2c u1 u3 1 v=C
B=q0 0 B=q0 0 0 1 0
0 0 0 0 0 0 1
2 66 66 66 66 66 64
3 77 77 77 77 77 75
;
and for matrix A3,
KA3¼ diag k ð3Þ1 ;kð3Þ2 ; . . . ;kð3Þ7
¼ diagðu3 c; u3; u3þ c; u3; . . . ; u3Þ;
RA3¼ r1ð3Þ; 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=C Hþ u3c u1 u2 1 v=C
B=q0 0 B=q0 0 0 1 0
0 0 0 0 0 0 1
2 66 66 66 66 66 64
3 77 77 77 77 77 75
;
AjrkðjÞ¼ kðjÞk rðjÞk , j¼ 1; 2; 3, and k ¼ 1; 2; . . . ; 7. Here c ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi cðp þ BÞ=q
p is the speed of sound of the fluid, and the other notations appeared in the above formulae are set by C¼ c 1, K ¼ CP3
j¼1u2j=2, H ¼ E þ p=q, and v¼ C½ðp þ c1B1Þ=ðc1 1Þ ðp þ c2B2Þ=ðc2 1Þ.
For the ease of the latter reference, it is customary to write (11) into a more compact expression by oq
otþXNd
j¼1
fj o oxj
; q
¼ 0; ð13Þ
where fj is taken as the vector-value function of the form
fj¼ o oxj
quj
; o oxj
qu1uj
þ pd1j
; . . . ; o oxj
quNduj
þ pdNdj
; o oxj
qEuj
þ puj
; o oxj
B q0quj
; ujoY oxj
T
; for j¼ 1; 2; . . . ; Nd. It is easy to see that the functions fjdefined above reduce to the standard flux functions for a single component flow.
To end this section, it is worth mentioning that if for some reasons the variables Y1q1 and Y2q2 are the preferable ones to be used in the algorithm, we may reformulate (11) into an analogous form of the five- equation model advocated by Allaire et al. [3] for compressible multicomponent flow problems. Since this 2-density formulation of the model system will not be considered in the present work, to avoid any possible confusions that may cause, we omit the full description of that model here.
4. Approximate Riemann solvers
Before describing numerical methods to solve (11), we pause to discuss the construction of the Riemann problem solutions in one space dimension which is one of the major steps in our barotropic two-fluid algorithm. If we consider the case Nd¼ 3 as an example, the problem to be solved now is to find the so- lution of (11) in the direction normal to one of the x1x2, x2x3, and x1x3planes, with piecewise constant data qL and qR to 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 qR well enough so that the solution of the Riemann problem would consist of genuinely nonlinear waves such as shock and rarefaction, and linearly degenerate wave such as contact discontinuity; there is no vacuum region occurring in the solution. To simplify the notation in the following discussion, rather than usingðu1; u2; u3Þ for the velocity field as before, we use a simpler version ðu; v; wÞ instead. Without loss of generality, we only look at the Riemann problem in the direction normal to the x2x3-plane.
4.1. Shock-only solver
To begin, we are interested in a popular shock-only (or called two-shock) approximation of the Riemann solver that ignores the possibility of rarefaction waves and simply construct a solution in which each pair of the states is connected along the Hugoniot locus for a shock (cf. [7,8]). In this approach, the key step is to find the midstate ðum; pmÞ in the u–p phase plane so that it can connect to ðuL; pLÞ by a 1-shock, and to ðuR; pRÞ by a 3-shock. It is well known that this is equivalent to solving the following nonlinear equation in an iterative manner for the pressure pm:
hðpmÞ ¼ umRðpmÞ umLðpmÞ ¼ 0: ð14Þ
Here umLand umRare the velocities defined by connecting the states along the 1-shock and 3-shock curves, respectively,
umLðpÞ ¼ uLp pL
MLðpÞ; umRðpÞ ¼ uRþp pR
MRðpÞ; ð15Þ
with Midenoting the Lagrangian shock speed, for i¼ L or R. In the current application of the pressure law (5), we may compute Miquite easily by evaluating the formula
Mi2ð Þ ¼p
ppi
q1i q1miðpÞ if Y ¼ 1 or Y ¼ 0 Ci2 1þ c2ciþ1
i
pþBi piþBi 1
h i
if Y 2 ð0; 1Þ;
8<
: ð16Þ
where Ci¼ qici is the Lagrangian sound speed, and qmi is the midstate density on the i side,
qmiðpÞ ¼ q0i ppþBi
0iþBi
1=ci
if Y ¼ 1 or Y ¼ 0;
q1i Mpp2 i iðpÞ
h i1
if Y 2 ð0; 1Þ:
8>
<
>: ð17Þ
Note that (16) and (17) are as a result derived from the Rankine–Hugoniot jump conditions across the shock waves (cf. [9]).
When applying a standard root-finding approach such as the secant method to (14), we have a 2-step iteration scheme as follows:
pðnþ1Þm ¼ pðnÞm pmðnÞ pðn1Þm uðnÞmL uðn1ÞmL
þ u ðnÞmR uðn1ÞmR
uðnÞmR uðnÞmL
h i
; ð18Þ
where uðnÞmi ¼ umi½pðnÞm, for i ¼ L or R, and n ¼ 1; 2; . . . (until convergence). With a suitable choice of the starting values pmð0Þ and pð1Þm, method (18) typically converges to the exact solution pm at a superlinear rate [22]. For gas dynamics, it is a common practice to set pmð0Þ and pð1Þm by
pð0Þm ¼pRCLþ pLCR ðuR uLÞCLCR CLþ CR
;
pð1Þm ¼pRMLð0Þþ pLMRð0Þ ðuR uLÞMLð0ÞMRð0Þ MLð0Þþ MRð0Þ ;
ð19Þ
where Mið0Þ¼ Mi½pmð0Þ. Having that, we may assign uð0ÞmL and uð0ÞmR by uð0ÞmL ¼ uLpð0Þm pL
CL
; uð0ÞmR¼ uRþpð0Þm pR
CR
;
and uð1ÞmL and uð1ÞmR according to (15). After a satisfactory convergence of the scheme, um can then be cal- culated based on the formula:
um¼pL pRþ uLMLðpmÞ þ uRMRðpmÞ MLðpmÞ þ MRðpmÞ :
We note that alternatively we may use a 1-step Newton method for the solution of (14). In this case, the iteration scheme reads as follows:
pðnþ1Þm ¼ pðnÞm ZLðnÞZRðnÞ ZLðnÞþ ZRðnÞhuðnÞmR
uðnÞmLi
; ð20Þ
where ZiðnÞ¼ Zi½pmðnÞ, for i ¼ L or R, and n ¼ 1; 2; . . . (until convergence). Again with (5) it is an easy matter to evaluate Zi by
Zi2ð Þ ¼p dpm
dum
¼
2C2miðpÞ
Mi2þCmi2ðpÞMi if Y ¼ 1 or Y ¼ 0;
2Mi2
Mi2þCi2Mi if Y 2 ð0; 1Þ;
8<
: ð21Þ
where Cmi¼ qmicmi. Here the initial guess of the iteration pmð0Þcan be taken in the same way as in (19) for the secant method. It is known that the rate of convergence of method (20) is quadratic, when pmð0Þis sufficiently close to pm(cf. [22]).
Fig. 1 shows a typical solution structure of the two-fluid Riemann problem considered here. Clearly, in a shock-only approximate solver, we replace the leftward-going rarefaction wave by a 1-shock and so the solution consists of three discontinuities moving at constant speeds. Here the propagation speed of each discontinuity is determined by
k1¼ umMLðpmÞ
qmLðpmÞ; k2¼ um; k3¼ umþMRðpmÞ
qmRðpmÞ; ð22Þ
with the jumps across each of them computed by the difference between the states to the left and right of the discontinuity,
W1¼ qmL qL; W2¼ qmR qmL; W3¼ qR qmR; ð23Þ
where qmiis calculated from (12b) using the data: qmi, um, pm, vi, wi, Bi, q0i, p0i, and Yi, for i¼ L or R. As usual, wave propagation methods (to be described in Section 5) are based on using these propagating discontinuities to update the cell averages in the cells neighboring each interface.
4.2. HLL-type solver
We are next concerned with a simple variant of the Riemann solver based on the work of Harten et al.
[16] for general hyperbolic systems of conservation laws, oq
otþ o
oxfðqÞ ¼ 0; ð24Þ
where q2 Rn is the vector of conserved variables for a system of n equations, and f is 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 kL and kRto the left and right, separating three constant states in the x–t space. If we assume further that kLand kRare known a priori by some simple estimates based on the local information of the wave speeds (cf. [10,52]), then by using the integral form of the conservation laws over a sufficiently large control volume½M; M ½0; T , for some M and T 2 R, it is
Fig. 1. (a) A typical solution structure of the Riemann problem. (b) Solution structure of the shock-only and HLLC approximate Riemann solvers. Note that in the original HLL solver the wave family associated with k2is not present in the solution.
an easy matter to find the constant state in the middle region, denoted by qm, as the average of the exact solution over the interval½T kL; TkR at time T ,
qm¼ 1
TðkR kLÞ Z TkR
TkL
qðx; T Þ dx ¼kRqR kLqL f ðqRÞ þ f ðqLÞ kR kL
; ð25Þ
where fðqiÞ is the flux evaluated at the state qi, for i¼ L or R.
While the above 2-wave HLL solver has been used quite successfully in many numerical methods for computing approximate solutions of (24), it is known that the numerical result obtained by using this solver is too diffusive for problems with contact discontinuities (cf. [52]). In addition to that, because of the lack of information on the structure of the interfaces, it is not an adequate approach at all for general multi- component problems.
To improve upon this 2-wave Riemann solver, we take a method suggested by Toro et al. [6,53] in that an additional middle wave of speed um is included 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 components of qm computed from (25), we may simply set um¼ qð2Þm=qð1Þm, where qðiÞm is the ith com- ponent of the vector qm, see [52] for the various other ways to compute um.
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 using the integral form of the conservation 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 applied to the current two-fluid model (11), it is not difficult to show that the result is
qð1Þmi ¼ fmið1Þ ki um
; qð2Þmi ¼ umqð1Þmi; qð3Þmi ¼ viqð1Þmi; qð4Þmi ¼ wiqð1Þmi; qð5Þmi ¼fmið5Þþ umumfmið1Þ fmið2Þ
ki um
; qð6Þmi ¼ fmið6Þ ki um
; qð7Þmi ¼ qð7Þi ;
ð26Þ
where fmiðiÞ¼ kiqðiÞi fðiÞðqiÞ represents the ith component of the vector fmi, for i¼ L or R. It is easy to check that qðiÞmL and qðiÞmR satisfy the basic consistency condition of the integral form of the conservation laws,
um kL
kR kL
qðiÞmLþ kR um
kR kL
qðiÞmR¼ qðiÞm;
for i¼ 1; 2; . . . ; 6. As in the shock-only solver, we then set the speed of the three moving discontinuities by k1¼ kL, k2¼ um, k3¼ kR, with the choice of kL ¼ minðuL cL; uR cRÞ and kR ¼ maxðuLþ cL; uRþ cRÞ as an example, and the jumps across each of them by (23) using the result (26).
To end this section, it should be mentioned that, when computing the Riemann problem solution of (11), the iterative shock-only solver described in the beginning of this section is a much more expensive method to use as compared to the non-iterative 3-wave HLLC solver. But due to the careful construction of the middle state solution for the contact discontinuity, which is extremely important for many difficult multicomponent problems with strong shock waves and stiff equations of state, we observe typically a better behavior of the results when the shock-only solver is in use than the HLLC solver is in use in a numerical method for ap- proximating solutions of practical problems. Surely, to improve efficiency, one may take a hybrid approach that utilizes the HLLC solver only in case for a single-fluid Riemann problem, but employs the shock-only solver in case for two-fluid Riemann problems otherwise. Finally, we note that it is still unclear on how to define the associated average states in a suitable way that linearizes the matrices in (12a), leading to an efficient and accurate Roe solver for the class of two-fluid Riemann problem considered here.
5. Numerical methods
To find approximate solutions of our model system (11) for barotropic two-fluid problems, we use a high-resolution wave propagation method developed by LeVeque [25,27] for general hyperbolic systems of partial differential equations. This method is a variant of the fluctuation-and-signal scheme of Roe [39,40]
in that we solve one-dimensional Riemann problems at each cell interface, and use the resulting waves (i.e., discontinuities moving at constant speeds) to update the solutions in neighboring grid cells. To achieve high resolution (i.e., second-order accurate on smooth solutions, and sharp and monotone profiles on discon- tinuous solutions), we introduce slopes and limiters to the method as in many other high-resolution schemes for conservation laws [14,26,51]. Here, for simplicity, rather than using x1, x2, and x3for the spatial vari- ables, we take the often-employed x, y, and z in a respective manner, instead.
5.1. One-dimensional case
We begin our discussion by reviewing the basic idea of the method for solving the simplest Nd¼ 1 case of our model system,
oq otþ f1
o ox; q
¼ 0;
with q and f1defined as in (13). To discretize the solution state, we take a uniform grid with fixed mesh spacing Dx in the x-direction and use a standard finite-volume formulation in which the value Qnj ap- proximates the cell average of the solution over the grid cell½xj; xjþ1 at time tn,
Qnj 1 Dx
Z xjþ1 xj
qðx; tnÞ dx:
The time step from the current time tn to the next tnþ1 is denoted by Dt.
5.1.1. First-order method
In this numerical discretization setup, a first-order accurate version of the method in wave-propagation form is a Godunov-type scheme that can be written as
Qnþ1j ¼ QnjDt Dx
Xmw
m¼1
kmWm
n
jþ1þ kþmWm
n
j; ð27Þ
where km and Wm are solutions of the mth wave family, for m¼ 1; 2; . . . ; mw, obtained from solving the Riemann problems at cell interfaces xjand xjþ1with the properly chosen approximate solver as described in Section 4, and k ¼ minðk; 0Þ, kþ¼ maxðk; 0Þ. It is known that method (27) belongs to a class of upwind schemes and is stable when the typical CFL (Courant–Friedrichs–Lewy) condition is satisfied (cf.
[14,26,28]). Moreover, it is not difficult to show that the method is quasi-conservative in the sense that when applying the method to (11) not only the conservation laws but also the transport equations are approx- imated in a consistent manner by the method, see [46] for the details.
In devising an efficient solver for multicomponent problems, it is of fundamental importance to look into an interface-only problem (cf. [45,46]) and check to see whether or not an oscillation-free results can be obtained by using the proposed numerical algorithm. Without loss of generality, here we consider a single two-fluid Riemann problem where at cell interface xj the initial data consists of uniform pressure p0 and constant velocity u0to the left and right of the interface, but having jumps on the other state variables such as q, c, B, q0, and Y . If we now assume a positive velocity u0>0, and take a shock-only solver (see Section 4.1) for the Riemann problem solution, from (27) the cell average Qnj would be updated by
Qnþ1j ¼ QnjDt
Dxðk2W2Þnj; ð28aÞ
or equivalently by q qu qE Bq=q0
Y 2 66 66 4
3 77 77 5
nþ1
j
¼ q qu0
qE Bq=q0
Y 2 66 66 4
3 77 77 5
n
j
Dt Dxu0
Dq u0Dq DðqEÞ DðBq=q0Þ
DY 2 66 66 4
3 77 77 5
n
j
; ð28bÞ
when expressing (28a) in terms of the solution states of the problem. Noting that in this case the difference operator D is simply applied to the Riemann data Qnj1 and Qnj on the left and right of the interface.
If we substitute the first equation of (28b) into the second one, it follows quite easily that we have the expected quantity of the particle velocity unþ1j ¼ u0. With this result, the third equation of (28b) can be simplified to
ðqeÞnþ1j ¼ ðqeÞnjDt
Dxu0DðqeÞnj or alternatively to
pþ cB c 1
B q0q
nþ1 j
¼ pþ cB c 1
B q0q
n j
Dt
Dxu0D pþ cB c 1
B q0q
n j
;
when employing the equation of state of the form (4). It is important to note that, when the fluid is barotropic, the total internal energy of the flow can be calculated directly by using (4) together with (2) for the pressure term. Thus we may use it to define the total energy for purely barotropic flows as in the case for the non-barotropic flows (i.e., the mixture of two different barotropic fluid components), yielding the vi- ability of the above formula, no matter what type of the fluid component is present in the Riemann data.
If we proceed our computation by applying the fourth equation of (28b) to the above equation, we obtain further simplification
pþ cB c 1
nþ1
j
¼ pþ cB c 1
n
j
Dt
Dxu0D pþ cB c 1
n
j
:
Now the replacement of c and B by the volume-fraction relations in (9b) would lead to an alternative form of the above equation as
X2
i¼1
Yipþ ciBi
ci 1
!nþ1
j
¼ X2
i¼1
p0þ ciBi
ci 1 ð ÞYi nj
!
Dt
Dxu0 X2
i¼1
p0þ ciBi
ci 1 ðDYiÞnj
! :
Then, after a simple manipulation, we get the desired pressure equilibrium pnþ1j ¼ p0, under the condition of the fulfillment of the difference equations for the volume fraction Yi,
Yi
ð Þnþ1j ¼ Yð Þi njDt
Dxu0D Yð Þi nj
for i¼ 1; 2. Note that this is exactly the last equation of (28b), because we have employed the notations in which Y1¼ Y and Y2¼ 1 Y .
Having shown this, it is also easy to demonstrate that, if the HLLC solver (see Section 4.2) is used instead for solving the aforementioned interface-only problem, we get the same Riemann problem solution as presented in (28b), and hence obtain the same numerical result of the problem.
5.1.2. High-resolution method
To improve the accuracy of method (27) to second-order, it is a customary approach by first introducing correction waves in a piecewise-linear form with zero mean value and then propagating each wave over the time step Dt to update the cell averages it overlaps. Without providing the details here (cf. [28,29] for example), with the corrections, (27) is modified by
Qnþ1j :¼ Qnþ1j Dt 2Dx
Xmw
m¼1
jkmj 1
jkmjDt Dx
Wfm
n jþ1
jkmj 1
jkmjDt Dx
Wfm
n j
; ð29Þ
where fWm is a limited value of Wm obtained by comparing Wm with the corresponding Wm from the neighboring Riemann problem to the left (if km>0) or to the right (if km<0).
Note that with the use of the wave-form representation of the Riemann solver as described in Section 4, it is quite common to perform the limiting procedure over each component of the wave via a limiter function U (e.g., by using the minmod function UðhÞ ¼ maxð0; minð1; hÞÞ or some others as discussed in [51]), and set
WfðiÞmj¼ U h ðiÞmj
WðiÞmj with hðiÞmj¼WðiÞmJ
WðiÞmj; J¼ j 1 if kmjP0;
jþ 1 if kmj<0;
ð30Þ
where WðiÞmj is the ith component of Wmj. While this approach works in a satisfactory manner without causing any wrong oscillations for the 1- and 3-waves, it was mentioned in [46,48] that, for the 2-wave, the third limited component on the total energy fWð3Þ2j should be modified to ensure a consistent approximation of that term, yielding the desired pressure equilibrium near the interfaces.
To demonstrate how the modification needs to be done in the current barotropic flow model (11), we use the interface-only problem described in Section 5.1.1 as an example. Then from (29) the cell average Qnþ1j obtained using (28a) would be updated by
Qnþ1j :¼ Qnþ1j Dt 2Dxjk2j 1
jk2jDt Dx
Wf2;jþ1
fW2j
n
; ð31aÞ
or equivalently by
q qu qE Bq=q0
Y 2 66 66 66 4
3 77 77 77 5
nþ1
j
:¼ q qu qE Bq=q0
Y 2 66 66 66 4
3 77 77 77 5
nþ1
j
l 2ð lÞ1
UhðDqÞj=ðDqÞjþ1i ðDqÞjþ1 u0UhðDqÞj=ðDqÞjþ1i
ðDqÞjþ1 UhðDqEÞj=ðDqEÞjþ1i
ðDqEÞjþ1 UhðDBq=q0Þj=ðDBq=q0Þjþ1i
ðDBq=q0Þjþ1 UhðDY Þj=ðDY Þjþ1i
ðDY Þjþ1 2
66 66 66 66 66 66 4
3 77 77 77 77 77 77 5
n
þl 2ð lÞ1
UhðDqÞj1=ðDqÞji ðDqÞj u0UhðDqÞj1=ðDqÞji
ðDqÞj UhðDqEÞj1=ðDqEÞji
ðDqEÞj UhðDBq=q0Þj1=ðDBq=q0Þji
ðDBq=q0Þj UhðDY Þj1=ðDY Þji
ðDY Þj 2
66 66 66 66 66 66 4
3 77 77 77 77 77 77 5
n
; ð31bÞ