### A ﬂuid-mixture type algorithm for barotropic two-ﬂuid ﬂow 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-ﬂuid ﬂow 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 ﬂuid 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 diﬀerent barotropic ﬂuid 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 ﬂow; Tait equation of state; Stiﬀened gas equation of state; Fluid-mixture model; Wave propagation method

1. Introduction

In this paper, we are concerned with a model barotropic two-ﬂuid ﬂow problem where the ﬂow 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

þX^{N}^{d}

j¼1

o
ox_{j}

quj

quiu_{j}þ pðqÞdij

¼ 0 for i¼ 1; 2; . . . Nd ð1Þ

for the basic equations of motion of each ﬂuid 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 ﬂuid 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
q_{0}

^{c}

B: ð2Þ

Here p0and q_{0}are 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-ﬂuid 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 ﬂow problem numerically, we want to use a generalization of the classical shock-capturing method designed originally for single component ﬂows. For non-barotropic multicom- ponent problems, it is known in the literature that the principal diﬃculty in the usual extension is the occurrence of spurious pressure oscillations when two or more ﬂuid 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 ﬂows, 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 ﬁrst introduce a non-isentropic form of the Tait equation of state as a basis to the modeling of the mixing between two diﬀerent barotropic ﬂuid components. Then with the help of a volume-fraction function, we deﬁne a hybrid equation of state so that the pressure of the ﬂuid can be determined explicitly no matter what ﬂuid 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 fulﬁllment 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-ﬂuid ﬂows with a Tait equation of state. Extension of the algorithm to other barotropic ﬂow 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-ﬂow 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 diﬀerent barotropic ﬂows within a grid cell. In Section 3, we describe the model equations that is easy to use for numerical approximation of barotropic two-ﬂuid 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 ﬂuid-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 diﬀerent ﬂuid components within a grid cell. In the current instance of barotropic
ﬂows, the basic idea of the proposed method is simple. We view the mixture of two diﬀerent barotropic
ﬂuids, where each of them is described by the same Tait equation of state (2) but possibly with a diﬀerent set
of material-dependent quantities: c, B, and q_{0}, as a non-barotropic ﬂuid, and use an entropy-dependent
equation of state for describing the thermodynamic behavior of the ﬂuid mixture instead. By using the ﬁrst
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 speciﬁc entropy S as:

pðq; SÞ ¼ AðSÞ pð 0þ BÞ q
q_{0}

c

B; ð3Þ

where AðSÞ ¼ exp½ðS S0Þ=CV and CVdenotes the speciﬁc heat at constant volume. Clearly, (3) reduces to the barotropic ﬂow 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 diﬀerent barotropic ﬂuids. Numerical experiences indicate that when DS is not large, i.e., entropy variation of the ﬂuids 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

q_{0}

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 ﬂuid dynamical literature an equation of state of the form (4) is typically called the stiﬀened gas or Tammann equation of state (cf. [14,17]). It is also easy to show that the internal energy of the purely barotropic ﬂuid 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 ﬂuid component within a grid cell (this is a standard way to do in many two-phase ﬂow solver [42,56,58]).

For instance, when grid cells contain only ﬂuid-component 1 we may set Y ¼ 1, and so when grid cells contain only ﬂuid-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 ﬂuid-components 1 and 2 occupied by the volume fractions Y and 1 Y , respectively. With this deﬁnition of Y , the pressure of a barotropic two-ﬂuid ﬂow problem in all the ﬂuid- component scenarios within a grid cell can be determined straightforwardly by:

p¼

p01þ B1

ð Þ _{q}^{q}

01

c_{1}

B1 if Y ¼ 1 ðbarotropic fluid 1Þ;

p02þ B2

ð Þ _{q}^{q}

02

c_{2}

B2 if Y ¼ 0 ðbarotropic fluid 2Þ;

ðc 1Þq e þ^{B}_{q}

0

cB if Y 2 ð0; 1Þ ðmixture of two barotropic fluidsÞ 8>

>>

<

>>

>:

ð5Þ

provided that all the variables appeared there are deﬁned 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 ﬂuid belongs to a set of real numbers. Certainly, it is both interesting and important to include the cavitation and phase-transition eﬀects to the present constitutive model in a region where the pressure drops to a critical value, but this subject matter is a very diﬃcult 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-ﬂuid 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 ﬂow. In regions where Y 2 ð0; 1Þ, however, since the mixture of two barotropic ﬂuid components is modeled in a non-isentropic manner, see Section 2, the original isentropic Eq. (1) are not suﬃcient to this instance, and so some supplementary equations need to be considered to provide further information of the ﬂuid mixture.

It is apparent that, because the thermodynamic behavior of the ﬂuid mixture depends on the entropy as we have postulated, the conservation law for the total energy of the form

o

otðqEÞ þX^{N}^{d}

j¼1

o oxj

qEuj

þ puj

¼ 0 ð6Þ

should be incorporated in the model system, where E¼ e þPNd

j¼1u^{2}_{j}=2 is the speciﬁc 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 eﬀective equations for the problem-dependent material quantities
so that the pressure can be computed easily from the equation of state.

3.1. c-based eﬀective equations

To derive the aforementioned eﬀective equations for the mixture of material quantities in the present stiﬀened 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 þX^{N}^{d}

j¼1

u_{j}oq
oxj

¼ 0;

o

otðqeÞ þX^{N}^{d}

j¼1

u_{j} o
oxj

ðqeÞ ¼ 0;

in a respective manner. By inserting the equation of state (4) into the latter one, we ﬁnd an alternative description of the energy equation

o ot

pþ cB c 1

B
q_{0}q

þX^{N}^{d}

j¼1

u_{j} o
oxj

pþ cB c 1

B
q_{0}q

¼ 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 ﬂuid mixtures of 1=ðc 1Þ and
cB=ðc 1Þ Bq=q_{0}as

o ot

1 c 1

þX^{N}^{d}

j¼1

u_{j} o
oxj

1 c 1

¼ 0; ð8aÞ

o ot

cB c 1

B
q_{0}q

þX^{N}^{d}

j¼1

u_{j} o
oxj

cB c 1

B
q_{0}q

¼ 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 satisﬁed 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 q_{0}. In our approach (cf. [46,48]), this is done by simply splitting (8b) into the following two parts:

o ot

cB c 1

þX^{N}^{d}

j¼1

u_{j} o
oxj

cB c 1

¼ 0; ð8cÞ

o ot

B
q_{0}q

þX^{N}^{d}

j¼1

u_{j} o
oxj

B
q_{0}q

¼ 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 ﬂuids near the interface. With that, there is no diﬃculty to compute the pressure from (4) as

p¼ qE

"

PNd

j¼1ðqujÞ^{2}

2q þB

q_{0}q 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 q_{0} 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 ﬂuid mixture the primitive form of (8d) should be modiﬁed by

o ot

B
q_{0}q

þX^{N}^{d}

j¼1

o oxj

B
q_{0}quj

¼ 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 eﬀective

equations for the mixture of the material quantities of the stiﬀened gas equation of state to be distinct from the other one presented below.

3.2. Y -based eﬀective equations

Before proceeding further, we note that to ﬁnd the initial ﬂuid mixtures, 1=ðc 1Þ, cB=ðc 1Þ, and Bq=q_{0},
that is necessary when we initialize the data for multicomponent ﬂow 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

q_{0}q¼ qe ¼X^{2}

i¼1

Y_{i}q_{i}e_{i}¼X^{2}

i¼1

Y_{i} p_{i}þ ciBi

c_{i} 1

Bi

q_{0i}q_{i}

: ð9aÞ

Here the subscript ‘‘i’’ denotes the state variable of ﬂuid component i. By taking a similar approach as employed in Section 3.1 for the derivation of the c-based eﬀective equation (cf. [45,46] also), it comes out easily a splitting of (9a) into the following set of relations:

1

c 1¼X^{2}

i¼1

Y_{i}

c_{i} 1; cB
c 1¼X^{2}

i¼1

Yi

c_{i}Bi

c_{i} 1; and B

q_{0}q¼X^{2}

i¼1

Yi

Bi

q_{0i}q_{i}; ð9bÞ

where in the process of splitting the terms the pressure p is chosen to satisfy the relation as p

c 1¼X^{2}

i¼1

Y_{i} p_{i}

c_{i} 1: ð9cÞ

With the ﬁrst 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 diﬀerent way to derive the same result.

Now with the volume-fraction notion of the states 1=ðc 1Þ, cB=ðc 1Þ, and Bq=q0being deﬁned by (9b), the c-based eﬀective equations may be rewritten into the form:

o ot

X^{2}

i¼1

Y_{i}
c_{i} 1

!
þX^{N}^{d}

j¼1

uj

o
ox_{j}

X^{2}

i¼1

Y_{i}
c_{i} 1

!

¼ 0; ð10aÞ

o ot

X^{2}

i¼1

Y_{i} c_{i}Bi

c_{i} 1

!
þX^{N}^{d}

j¼1

u_{j} o
oxj

X^{2}

i¼1

Y_{i} c_{i}Bi

c_{i} 1

!

¼ 0; ð10bÞ

o ot

X^{2}

i¼1

Y_{i}Bi

q_{0i}q_{i}

!
þX^{N}^{d}

j¼1

o
ox_{j}

X^{2}

i¼1

Y_{i}Bi

q_{0i}q_{i}u_{j}

!

¼ 0: ð10cÞ

After some simple algebraic manipulations (cf. [45]), from both (10a) and (10b), we ﬁnd the transport equation for each volume fraction Yi as

oYi

ot þX^{N}^{d}

j¼1

u_{j}oYi

oxj

¼ 0; ð10dÞ

while from (10c) we ﬁnd the mass conservation equation for each ﬂuid component i as o

otðY_{i}q_{i}Þ þX^{N}^{d}

j¼1

o oxj

Y_{i}q_{i}u_{j}

¼ 0; ð10eÞ

for i¼ 1; 2. Clearly, when the set of Yiand Yiq_{i}are 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 eﬀective 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 ﬂuid 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 aﬀecting anything, to a single one for Y as

oY
ot þX^{N}^{d}

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 ﬂuid-mixture type algorithm for multicomponent problems, one common practice (cf. [45,46,48])
is to consider the mixture of the total density, q¼ Y1q_{1}þ Y2q_{2}, as one of the basic variables, but is not to
use the mixtures of the separate ﬂuid densities, Y1q_{1}and Y2q_{2}, 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-ﬂuid ﬂow problems in more than one space dimension take the form:

oq
otþX^{N}^{d}

j¼1

o oxj

ðqujÞ ¼ 0;

o

otðquiÞ þX^{N}^{d}

j¼1

o oxj

ðquiu_{j}þ pdijÞ ¼ 0 for i¼ 1; 2; . . . ; Nd;
o

otðqEÞ þX^{N}^{d}

j¼1

o oxj

qEuj

þ puj

¼ 0 if Y 2 ð0; 1Þ;

o ot

B
q_{0}q

þX^{N}^{d}

j¼1

o oxj

B
q_{0}quj

¼ 0 if Y 2 ð0; 1Þ;

oY
ot þX^{N}^{d}

j¼1

u_{j}oY
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 ﬁrst 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 ﬂuid-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 ﬁnd c and B according to the ﬁrst two re- lations in (9b)). On the other instance, when Y ¼ 0 or Y ¼ 1, we just have the ﬁrst Ndþ 1 equations

governing the single-component barotropic ﬂow 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 ﬂuid-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 ﬂow are all in the region of the thermodynamic stability (this is the case we are interested in here), it is not diﬃcult 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þX^{N}^{d}

j¼1

AjðqÞoq

ox_{j}¼ 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) deﬁned by

q¼ q; qu1;qu2;qu3;qE;B
q_{0}q; Y

T

; ð12bÞ

and the matrices Aj, for j¼ 1; 2; 3, deﬁned by

A_{1}¼

0 1 0 0 0 0 0

K u^{2}_{1} u_{1}ð2 CÞ u2C u3C C C v

u1u_{2} u_{2} u_{1} 0 0 0 0

u1u_{3} u_{3} 0 u_{1} 0 0 0

u1ðK H Þ H u^{2}_{1}C u1u2C u1u3C u1ðC þ 1Þ u1C u1v

u1B=q_{0} B=q_{0} 0 0 0 u_{1} 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 u^{2}_{2} u1C u_{2}ð2 CÞ u3C C C v

u2u3 0 u3 u2 0 0 0

u_{2}ðK H Þ u1u_{2}C H u^{2}_{2}C u2u_{3}C u_{2}ðC þ 1Þ u2C u_{2}v

u2B=q0 0 B=q0 0 0 u_{2} 0

0 0 0 0 0 0 u_{2}

2 66 66 66 66 66 64

3 77 77 77 77 77 75

;

A_{3}¼

0 0 0 1 0 0 0

u1u_{3} u_{3} 0 u_{1} 0 0 0

u2u_{3} 0 u_{3} u_{2} 0 0 0

K u^{2}_{3} u1C u2C u_{3}ð2 CÞ C C v
u_{3}ðK H Þ u1u_{3}C u2u_{3}C H u^{2}_{3}C u_{3}ðC þ 1Þ u3C u_{3}v

u3B=q_{0} 0 0 B=q_{0} 0 u3 0

0 0 0 0 0 0 u_{3}

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; u_{1}þ c; u1; . . . ; u_{1}Þ;

R_{A}_{1}¼ 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

u_{2} u_{2} u_{2} 1 0 0 0

u_{3} u_{3} u_{3} 0 1 0 0

H u1c K=C Hþ u1c u_{2} u_{3} 1 v=C

B=q_{0} 0 B=q_{0} 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; u_{2}þ c; u2; . . . ; u_{2}Þ;

R_{A}_{1}¼ r_{1}^{ð2Þ}; r^{ð2Þ}_{2} ; . . . ; r^{ð2Þ}_{7}

¼

1 1 1 0 0 0 0

u1 u1 u1 1 0 0 0

u_{2} c u_{2} u_{2}þ c 0 0 0 0

u3 u3 u3 0 1 0 0

H u2c K=C Hþ u2c u_{1} u_{3} 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; u_{3}þ c; u3; . . . ; u_{3}Þ;

R_{A}_{3}¼ r_{1}^{ð3Þ}; r^{ð3Þ}_{2} ; . . . ; r^{ð3Þ}_{7}

¼

1 1 1 0 0 0 0

u1 u1 u1 1 0 0 0

u_{2} u_{2} u_{2} 0 1 0 0

u3 c u3 u3þ c 0 0 0 0
H u3c K=C Hþ u3c u_{1} u_{2} 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

;

A_{j}r_{k}^{ðjÞ}¼ k^{ðjÞ}_{k} r^{ðjÞ}_{k} , j¼ 1; 2; 3, and k ¼ 1; 2; . . . ; 7. Here c ¼ ﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
cðp þ BÞ=q

p is the speed of sound of the ﬂuid, and the other notations appeared in the above formulae are set by C¼ c 1, K ¼ CP3

j¼1u^{2}_{j}=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þX^{N}^{d}

j¼1

f_{j} o
oxj

; q

¼ 0; ð13Þ

where fj is taken as the vector-value function of the form

f_{j}¼ o
oxj

quj

; o oxj

qu1u_{j}

þ pd1j

; . . . ; o oxj

quNdu_{j}

þ pdNdj

; o oxj

qEuj

þ puj

; o oxj

B
q_{0}quj

; u_{j}oY
oxj

T

; for j¼ 1; 2; . . . ; Nd. It is easy to see that the functions fjdeﬁned above reduce to the standard ﬂux functions for a single component ﬂow.

To end this section, it is worth mentioning that if for some reasons the variables Y1q_{1} and Y2q_{2} are the
preferable ones to be used in the algorithm, we may reformulate (11) into an analogous form of the ﬁve-
equation model advocated by Allaire et al. [3] for compressible multicomponent ﬂow 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-ﬂuid
algorithm. If we consider the case Nd¼ 3 as an example, the problem to be solved now is to ﬁnd the so-
lution of (11) in the direction normal to one of the x1x_{2}, x2x_{3}, and x1x_{3}planes, with piecewise constant data
q_{L} 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; u_{2}; u_{3}Þ for the velocity ﬁeld 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
ﬁnd the midstate ðum; p_{m}Þ in the u–p phase plane so that it can connect to ðuL; p_{L}Þ 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 deﬁned 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

M_{i}^{2}ð Þ ¼p

ppi

q^{1}_{i} q^{1}miðpÞ if Y ¼ 1 or Y ¼ 0
C_{i}^{2} 1þ ^{c}_{2c}^{i}^{þ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 q_{mi} is the midstate density on the i side,

q_{mi}ðpÞ ¼ q_{0i} _{p}^{pþB}^{i}

0iþBi

1=c_{i}

if Y ¼ 1 or Y ¼ 0;

q^{1}_{i} _{M}^{pp}2 ^{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-ﬁnding approach such as the secant method to (14), we have a 2-step iteration scheme as follows:

p^{ðnþ1Þ}_{m} ¼ p^{ðnÞ}_{m} p_{m}^{ð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 p_{m}^{ð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 p_{m}^{ð0Þ} and p^{ð1Þ}_{m} by

p^{ð0Þ}_{m} ¼p_{R}C_{L}þ pLC_{R} ðuR uLÞCLC_{R}
C_{L}þ CR

;

p^{ð1Þ}_{m} ¼p_{R}M_{L}^{ð0Þ}þ pLM_{R}^{ð0Þ} ðuR uLÞM_{L}^{ð0Þ}M_{R}^{ð0Þ}
M_{L}^{ð0Þ}þ M_{R}^{ð0Þ} ;

ð19Þ

where M_{i}^{ð0Þ}¼ Mi½p_{m}^{ð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:

u_{m}¼p_{L} pRþ uLM_{L}ðpmÞ þ uRM_{R}ðpmÞ
M_{L}ð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} Z_{L}^{ðnÞ}Z_{R}^{ðnÞ}
Z_{L}^{ðnÞ}þ Z_{R}^{ðnÞ}hu^{ðnÞ}_{mR}

u^{ðnÞ}_{mL}i

; ð20Þ

where Z_{i}^{ðnÞ}¼ Zi½p_{m}^{ðnÞ}, for i ¼ L or R, and n ¼ 1; 2; . . . (until convergence). Again with (5) it is an easy matter
to evaluate Zi by

Z_{i}^{2}ð Þ ¼p dpm

dum

¼

2C^{2}_{mi}ðpÞ

M_{i}^{2}þCmi^{2}ðpÞM_{i} if Y ¼ 1 or Y ¼ 0;

2M_{i}^{2}

M_{i}^{2}þCi^{2}M_{i} if Y 2 ð0; 1Þ;

8<

: ð21Þ

where Cmi¼ qmic_{mi}. Here the initial guess of the iteration p_{m}^{ð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 p_{m}^{ð0Þ}is suﬃciently
close to pm(cf. [22]).

Fig. 1 shows a typical solution structure of the two-ﬂuid 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¼ umM_{L}ðpmÞ

q_{mL}ðpmÞ; k2¼ um; k3¼ umþM_{R}ðpmÞ

q_{mR}ðpmÞ; ð22Þ

with the jumps across each of them computed by the diﬀerence 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: q_{mi}, um, pm, vi, wi, Bi, q_{0i}, 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 R^{n} is the vector of conserved variables for a system of n equations, and f is the ﬂux 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 suﬃciently 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 ﬁnd 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 Tk_{R}

Tk_{L}

qðx; T Þ dx ¼kRq_{R} kLq_{L} f ðqRÞ þ f ðqLÞ
kR kL

; ð25Þ

where fðqiÞ is the ﬂux 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 diﬀusive 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 ﬁrst
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 ﬁnd 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-ﬂuid model (11), it is not diﬃcult to show that the result is

q^{ð1Þ}_{mi} ¼ f_{mi}^{ð1Þ}
ki um

; q^{ð2Þ}_{mi} ¼ umq^{ð1Þ}_{mi}; q^{ð3Þ}_{mi} ¼ viq^{ð1Þ}_{mi}; q^{ð4Þ}_{mi} ¼ wiq^{ð1Þ}_{mi};
q^{ð5Þ}_{mi} ¼f_{mi}^{ð5Þ}þ umumf_{mi}^{ð1Þ} f_{mi}^{ð2Þ}

ki um

; q^{ð6Þ}_{mi} ¼ f_{mi}^{ð6Þ}
ki um

; q^{ð7Þ}_{mi} ¼ q^{ð7Þ}_{i} ;

ð26Þ

where f_{mi}^{ð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,

u_{m} 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; u_{R} cRÞ and kR ¼ maxðuLþ cL; u_{R}þ 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 diﬃcult multicomponent problems with strong shock waves and stiﬀ 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 eﬃciency, one may take a hybrid approach that utilizes the HLLC solver only in case for a single-ﬂuid Riemann problem, but employs the shock-only solver in case for two-ﬂuid Riemann problems otherwise. Finally, we note that it is still unclear on how to deﬁne the associated average states in a suitable way that linearizes the matrices in (12a), leading to an eﬃcient and accurate Roe solver for the class of two-ﬂuid Riemann problem considered here.

5. Numerical methods

To ﬁnd approximate solutions of our model system (11) for barotropic two-ﬂuid problems, we use a high-resolution wave propagation method developed by LeVeque [25,27] for general hyperbolic systems of partial diﬀerential equations. This method is a variant of the ﬂuctuation-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 proﬁles 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 f1deﬁned as in (13). To discretize the solution state, we take a uniform grid with ﬁxed mesh
spacing Dx in the x-direction and use a standard ﬁnite-volume formulation in which the value Q^{n}_{j} ap-
proximates the cell average of the solution over the grid cell½xj; x_{jþ1} at time tn,

Q^{n}_{j} 1
Dx

Z xjþ1 xj

qðx; tnÞ dx:

The time step from the current time tn to the next t_{nþ1} is denoted by Dt.

5.1.1. First-order method

In this numerical discretization setup, a ﬁrst-order accurate version of the method in wave-propagation form is a Godunov-type scheme that can be written as

Q^{nþ1}_{j} ¼ Q^{n}_{j}Dt
Dx

X^{m}^{w}

m¼1

k^{}_{m}Wm

n

jþ1þ k^{þ}_{m}Wm

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 x_{j}and x_{jþ1}with 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 satisﬁed (cf.

[14,26,28]). Moreover, it is not diﬃcult 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 eﬃcient 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-ﬂuid 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 Q^{n}_{j} would be updated by

Q^{nþ1}_{j} ¼ Q^{n}_{j}Dt

Dxðk2W2Þ^{n}_{j}; ð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 diﬀerence
operator D is simply applied to the Riemann data Q^{n}_{j1} and Q^{n}_{j} on the left and right of the interface.

If we substitute the ﬁrst equation of (28b) into the second one, it follows quite easily that we have the expected
quantity of the particle velocity u^{nþ1}_{j} ¼ u0. With this result, the third equation of (28b) can be simpliﬁed to

ðqeÞ^{nþ1}_{j} ¼ ðqeÞ^{n}_{j}Dt

Dxu0DðqeÞ^{n}_{j}
or alternatively to

pþ cB c 1

B
q_{0}q

nþ1 j

¼ pþ cB c 1

B
q_{0}q

n j

Dt

Dxu0D pþ cB c 1

B
q_{0}q

n j

;

when employing the equation of state of the form (4). It is important to note that, when the ﬂuid is barotropic, the total internal energy of the ﬂow can be calculated directly by using (4) together with (2) for the pressure term. Thus we may use it to deﬁne the total energy for purely barotropic ﬂows as in the case for the non-barotropic ﬂows (i.e., the mixture of two diﬀerent barotropic ﬂuid components), yielding the vi- ability of the above formula, no matter what type of the ﬂuid 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 simpliﬁcation

pþ cB c 1

nþ1

j

¼ pþ cB c 1

n

j

Dt

Dxu_{0}D 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

X^{2}

i¼1

Y_{i}pþ c_{i}Bi

c_{i} 1

!nþ1

j

¼ X^{2}

i¼1

p_{0}þ c_{i}Bi

c_{i} 1 ð ÞY_{i} ^{n}_{j}

!

Dt

Dxu_{0} X^{2}

i¼1

p_{0}þ c_{i}Bi

c_{i} 1 ðDYiÞ^{n}_{j}

! :

Then, after a simple manipulation, we get the desired pressure equilibrium p^{nþ1}_{j} ¼ p0, under the condition of
the fulﬁllment of the diﬀerence equations for the volume fraction Yi,

Y_{i}

ð Þ^{nþ1}_{j} ¼ Yð Þi ^{n}_{j}Dt

Dxu_{0}D Yð Þi ^{n}_{j}

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 ﬁrst 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 modiﬁed by

Q^{nþ1}_{j} :¼ Q^{nþ1}_{j} Dt
2Dx

X^{m}^{w}

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 k_{m}>0) or to the right (if k_{m}<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 modiﬁed to ensure a consistent approximation
of that term, yielding the desired pressure equilibrium near the interfaces.

To demonstrate how the modiﬁcation needs to be done in the current barotropic ﬂow model (11), we use
the interface-only problem described in Section 5.1.1 as an example. Then from (29) the cell average Q^{nþ1}_{j}
obtained using (28a) would be updated by

Q^{nþ1}_{j} :¼ Q^{nþ1}_{j} 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þ1}i
ðDqÞ_{jþ1}
u_{0}UhðDqÞ_{j}=ðDqÞ_{jþ1}i

ðDqÞ_{jþ1}
UhðDqEÞ_{j}=ðDqEÞ_{jþ1}i

ðDqEÞ_{jþ1}
UhðDBq=q_{0}Þ_{j}=ðDBq=q_{0}Þ_{jþ1}i

ðDBq=q_{0}Þ_{jþ1}
UhðDY Þ_{j}=ðDY Þ_{jþ1}i

ð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Þ_{j}i
ðDqÞ_{j}
u_{0}UhðDqÞ_{j1}=ðDqÞ_{j}i

ðDqÞ_{j}
UhðDqEÞ_{j1}=ðDqEÞ_{j}i

ðDqEÞ_{j}
UhðDBq=q0Þ_{j1}=ðDBq=q0Þ_{j}i

ðDBq=q0Þ_{j}
UhðDY Þ_{j1}=ðDY Þ_{j}i

ðDY Þ_{j}
2

66 66 66 66 66 66 4

3 77 77 77 77 77 77 5

n

; ð31bÞ