• 沒有找到結果。

Keh-MingShyue Afluid-mixturetypealgorithmforbarotropictwo-fluidflowproblems

N/A
N/A
Protected

Academic year: 2022

Share "Keh-MingShyue Afluid-mixturetypealgorithmforbarotropictwo-fluidflowproblems"

Copied!
31
0
0

加載中.... (立即查看全文)

全文

(1)

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

(2)

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.

(3)

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:

(4)

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;

(5)

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

(6)

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Þ

(7)

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

(8)

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 :

(9)

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

(10)

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 Ci2c2ciþ1

i

 

pþBi piþBi 1

 

h i

if Y 2 ð0; 1Þ;

8<

: ð16Þ

(11)

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]).

(12)

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.

(13)

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.

(14)

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

(15)

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.

(16)

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Þ

參考文獻

相關文件

• There are important problems for which there are no known efficient deterministic algorithms but for which very efficient randomized algorithms exist. – Extraction of square roots,

One of the main results is the bound on the vanishing order of a nontrivial solution to the Stokes system, which is a quantitative version of the strong unique continuation prop-

One of the main results is the bound on the vanishing order of a nontrivial solution u satisfying the Stokes system, which is a quantitative version of the strong unique

Reading Task 6: Genre Structure and Language Features. • Now let’s look at how language features (e.g. sentence patterns) are connected to the structure

• Examples of items NOT recognised for fee calculation*: staff gathering/ welfare/ meal allowances, expenses related to event celebrations without student participation,

Problems : Strong coupling, many body, solitons, …. Note: no need for theories to be

For problems 1 to 9 find the general solution and/or the particular solution that satisfy the given initial conditions:. For problems 11 to 14 find the order of the ODE and

Since it is so, what do we cultivate for?People are looking for the ways to improve the mental state, and the courage or wisdom to face the hard moments.. But the ways of improving