行政院國家科學委員會專題研究計畫 成果報告
可壓縮性多相位流與非凸性狀態方程的數值解法(2/2)
計畫類別： 個別型計畫
計畫編號： NSC92-2115-M-002-018-
執行期間： 92 年 08 月 01 日至 93 年 12 月 31 日
執行單位： 國立臺灣大學數學系暨研究所
計畫主持人： 薛克民
報告類型： 完整報告
報告附件： 出席國際會議研究心得報告及發表論文
處理方式： 本計畫可公開查詢
中 華 民 國 94 年 4 月 19 日
(will be inserted by the editor)
A volume-fraction based algorithm for hybrid barotropic and
non-barotropic two-fluid flow problems
Keh-Ming Shyue
Department of Mathematics, National Taiwan University, Taipei 106, Taiwan, Republic of China
Received
Abstract. The aim of this paper is to describe a simple Eulerian interface-capturing approach for the effi-cient numerical resolution of a hybrid barotropic and non-barotropic two-fluid flow problem in more than one space dimension. We use the compressible Euler equations as a model system with the thermodynamic property of each of the barotropic and non-barotropic fluid components characterized by the Tait and Noble-Abel equations of state, respectively. The algorithm is based on a volume fraction formulation of the equations together with an extended equation of state that is devised to give an approximate treatment for the mixture of more than one fluid component within a grid cell. A standard high-resolution wave propa-gation method is employed to solve the proposed two-fluid model with the dimensional-splitting technique incorporated in the method for multidimensional problems. Several numerical results are presented in one and two space dimensions that show the feasibility of the algorithm as applied to a reasonable class of practical problems without the occurrence of any spurious oscillation in the pressure near the smeared material interfaces. This includes, in particular, solutions for a study on the variation of the jet velocity with the incident shock pressure arising from the collapse of an air cavity in water under a shock wave. Key words: Volume-fraction based algorithm, Tait equation of state, Noble-Abel equation of state, Gas-cavity collapse, Jetting
PACS: 02.70.Bf, 46.15.-x, 47.11.+j, 47.55.Kf
1 Introduction
This paper is concerned with a simplified compressible two-phase flow problem where there is a material inter-face separating regions of two different barotropic and non-barotropic fluid components within an Nd≥ 1 spatial
domain. In this problem, the flow regime of interest is as-sumed to be homogeneous with no jumps in the pressure and velocity (the normal component of it) across the in-terface, and the physical effects such as the viscosity and heat conduction are assumed to be small, and hence can be ignored in the problem formulation. With that, we use an Eulerian viewpoint of the model system, that on the barotropic part of the domain, the fluid is governed by the isentropic version of the compressible Euler equations as
∂ ∂t ρ ρui + Nd X j=1 ∂ ∂xj ρuj ρuiuj+ p(ρ)δij = 0, (1)
for i = 1, 2, · · · , Nd, while on the non-barotropic part of
the domain, the fluid is governed by the full set of the
Correspondence to: K.-M. Shyue (Tel: 8862-3366-2866, Fax: 8862-2391-4439, e-mail: shyue@math.ntu.edu.tw)
compressible Euler equations as ∂ ∂t ρuρi ρE + Nd X j=1 ∂ ∂xj ρuiuρuj+ pδj ij ρEuj+ puj = 0. (2)
Here ρ, uj, p, E, and δij denote the density, the particle
velocity in the xj-direction, the pressure, the specific total
energy, and the Kronecker delta, respectively.
To complete the model, in this work, the constitutive law for each of the barotropic and non-barotropic com-ponents 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 [23]),
p(ρ) = (p0+ B)
ρ ρ0
γ
− B, (3) and the Noble-Abel equation of state for real gases (or called the constant covolume equation of state [37]),
p(ρ, e) = γ − 1 1 − bρ ρe, (4) in a respective manner. Note that in the above expres-sions e represents the internal energy, ρ0 and p0
pressure-like constant, b represents a finite size of volume occupied by molecules (0 ≤ b < 1/ρ), and γ is the adi-abatic constant (γ > 1), see [19,22,24,42] for a typical set of data of practical importance. As usual, we have E = e + PNd
j=1u2j/2. Sample applications of this
two-fluid model are, for instance, to the simulation of a shock-induced collapse of gas cavities in liquid [6,7], or to the simulation of an underwater explosion bubble [20,41].
For an efficient numerical resolution of this hybrid barotropic and non-barotropic two-fluid flow problem, we want to use a generalization of the classical shock-capturing method designed originally for single compo-nent flows in an Eulerian framework. It is known in the literature that the principal difficulty in the usual modi-fication is the occurrence of spurious pressure oscillations when two or more fluid components are present in a uniform Cartesian grid cell. In the case of a barotropic or non-barotropic multicomponent problem, there has been quite a few numerical methods developed for that matter, see [1,2,8,13,14,28,40] and the references cited therein for more exposition. In the current interest of the problems consisting of separate barotropic and non-barotropic regions, there is, however, a relatively few attempts devoted to the subject, see [4,20,34] for an example using a somewhat complicated Lagrangian-type approach as compared to the current Eulerian one.
The method we take here is an extension of the previ-ous work advocated by the author [30,31,33] in that a mix-ture equation of state is introduced first as a basis to the modeling of the mixing between two different barotropic and non-barotropic fluid components. Then with the help of a volume-fraction function, we define an extended equa-tion of state so that the pressure of the fluid can be deter-mined explicitly no matter what fluid component (pure or not) is within a grid cell (see Section 2.1). Having that, as before, we are able to derive a volume-fraction based of the model system that consists of the full Euler equations for the basic conserved variables and an additional set of equations for the volume fraction as well as the partial density of any one of the two fluid components in the prob-lem. In our approach, the latter equations are included in the algorithm primarily for an easy determination of the approximate location of the interface and also the com-putation of the pressure from the equation of state. It is important to note that these equations are put in a form so as to ensure a consistent modeling of the energy equa-tion near the smeared interfaces, and also the fulfillment of the mass equation in the other single component regions (see Section 2.2).
To find approximate solutions of our model system, we use a high-order Godunov method based on a wave-propagation formulation with the dimensional-splitting technique incorporated in the method for multidimen-sional problems. Numerical results to be presented in Section 4 show that this is a viable approach to a rea-sonable class of practical problems in one and two space dimensions, without producing any wrong oscillations in the pressure near the interfaces.
It should be emphasized that the methodology we have proposed here is by no means limited to the two-fluid flows with the above chosen equations of state. Extension of the algorithm to problems involving more than two fluid com-ponents and more complicated equations of state, can be made in a straightforward manner by following the idea described in this paper and the five equation model of Al-laire et al. [2]. Without going into the details for that, our goal is to establish the basic solution strategy and validate its use via some sample numerical experimentations.
The format of this paper is as follows. In Section 2, we describe the basic mathematical equations in a fluid-mixture form that will be used later in a numerical method for constructing approximate solutions of our model hy-brid barotropic and non-barotropic two-fluid problems. In Section 3, we give a brief review of the HLLC approximate Riemann solver and the high-resolution wave-propagation methods for solving problems in one and more than one space dimensions. Numerical results of some sample ex-amples are shown in Section 4.
2 Mathematical formulation
2.1 Equations of state
To start out, we want to build a mixture equation of state that is prerequisite in our approach for modeling the nu-merical mixing between the barotropic and non-barotropic fluid components within a grid cell. The first step, we take here, is to assume that thermodynamically the mixture of the fluids characterized by (3) and (4) behaves like a typ-ical non-barotropic fluid. Then using the first and second laws of thermodynamics, we may rewrite (3) and (4) in terms of the specific entropy S and the specific volume V = 1/ρ as p(V, S) = A(S) (p0+ B) V0 V γ − B, (5) and p(V, S) = A(S)p0 V0− b V − b γ , (6) for each in turn, where A(S) = R exp [(S − S0)/CV]
(cf. [9]). Here CV represents the specific heat at constant
volume, R is the universal gas constant, and S0 is the
specific entropy at a reference state. As in the previous work for multicomponent problems with a van der Waals equation of state [30], we assume further that all the fluid components under concerned is in an adiabatic equilib-rium with the same entropy, and accordingly we propose to define an equation of state for the fluid mixture that combines (5) with (6) to be the form
p(V, S) = A(S) (p0+ B) V0− b V − b γ − B. (7) Notice that, with the use of this so to be called modi-fied Noble-Abel equation of state (7), we recover not only
the case (6) when B → 0, but also the case (5) and so (3) as we set S = S0in (5), when b → 0. In the non-limiting case,
however, where none of the material-dependent quantities B and b is close to zero, (7) does give a way to the repre-sentation of the cases in between, that is to the mixing of the present barotropic and non-barotropic components.
It should be mentioned that, as a practical matter for many fluid flow problems away from cavitation (which is the situation studied here), rather than using (7) for the fundamental thermodynamic description of the fluid mixture, it is customary to formulate the equation in a way with the often-used variables in compressible hydro-dynamics: ρ and e instead. Again, by employing the basic thermodynamic principles, this can be done quite easily, leading to p(ρ, e) = γ − 1 1 − bρ (ρe − B) − B. (8) As to the computation of the fluid temperature T which is important in a wide variety of applications, we may use one of the formulas,
p(ρ, T ) = ρRT 1 − bρ− B, e(ρ, T ) = RT γ − 1+ B ρ, for realization.
Now, for an easy identification of the type of fluid com-ponent within a grid cell, a volume-fraction function Y is introduced here for that purpose. For example, when grid cells contain only the barotropic component we may set Y = 1, and so when grid cells contain only the non-barotropic component we set Y = 0. In case there are some cells cut by the interfaces where Y ∈ (0, 1), we then have both of these components occupied by the volume fractions Y and 1 − Y for each separately. With this defi-nition of Y , the pressure of our two-fluid flow problem in all the fluid-component scenarios within a grid cell can be determined straightforwardly by p = (p0+ B) ρ ρ0 γ − B if Y = 1 γ − 1 1 − bρ (ρe − B) − B if Y 6= 1, (9)
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 ma-terials of interest is limited by the stability requirement that the speed of sound of the fluid belongs to a set of real numbers.
2.2 Equations of motion
With the proposed constitutive law (9), it is clear that in regions where Y = 0 or Y = 1, there is no problem in
using (1) and (2), respectively, for the complete descrip-tion of the behavior of the underlying single-component flow. In regions where Y ∈ (0, 1), since the fluid mixture is modeled in a non-barotropic manner, it is quite com-mon to apply (2) as a model system that describes the motion of the mixtures of the basic conserved variables. Besides that, as in our previous work on numerical meth-ods for compressible multicomponent problems (cf. [29,30, 32,33]), we also bring in a set of effective equations for the problem-dependent material quantities so that the pres-sure can be computed easily from the equation of state.
2.2.1γ-based effective equations
To derive the aforementioned effective equations for the mixture of material quantities in the current hybrid barotropic and non-barotropic flow application, we start with an interface-only problem as usual where both the pressure and each component of the particle velocities are constants in the domain, while the other variables such as the density and the material quantities are having jumps across some interfaces. Then, from (2), it is easy to obtain equations for the time-dependent behavior of the density and total internal energy as
∂ρ ∂t + Nd X j=1 uj ∂ρ ∂xj = 0, ∂ ∂t(ρe) + Nd X j=1 uj ∂ ∂xj (ρe) = 0,
in a respective manner. Now, by inserting the modified Noble-Abel equation of state (8) into the latter one, we find an alternative description of the energy equation
∂ ∂t 1 − bρ γ − 1p + γ − bρ γ − 1 B + Nd X j=1 uj ∂ ∂xj 1 − bρ γ − 1p + γ − bρ γ − 1 B = 0 (10)
that is in relation to not only the pressure, but also the density and the material quantities: γ, b, and B.
In our algorithm, to maintain the pressure in equilib-rium as it should be for this interface-only problem, we split (10) into the following two equations for the fluid mixtures (1 − bρ)/(γ − 1) and (γ − bρ)B/(γ − 1) as ∂ ∂t 1 − bρ γ − 1 + Nd X j=1 uj ∂ ∂xj 1 − bρ γ − 1 = 0, (11) ∂ ∂t γ − bρ γ − 1B + Nd X j=1 uj ∂ ∂xj γ − bρ γ − 1B = 0, (12)
respectively. We emphasize that in order to have the cor-rect pressure equilibrium in (10), these are the two key
equations that should be satisfied and approximated con-sistently (when the problem is solved numerically, see Sec-tion 3). On the other hand, as a practical matter, it is obvious that, in addition to (11) and (12), we need to im-pose one more additional condition so as to have enough equations for the determination of all the three material quantities in (8). In our approach, see [30], this is done by simply breaking (11) into the following two parts:
∂ ∂t 1 γ − 1 + Nd X j=1 uj ∂ ∂xj 1 γ − 1 = 0, (13) ∂ ∂t bρ γ − 1 + Nd X j=1 uj ∂ ∂xj bρ γ − 1 = 0. (14)
Having done so, we arrive at a system of three equa-tions (12), (13), and (14) for the variables (γ−bρ)B/(γ−1), 1/(γ −1), and bρ/(γ −1), in the order mentioned. Combin-ing them to (2) yields a model system that is fundamen-tal in our algorithm for describing the behavior of the nu-merical mixing between the barotropic and non-barotropic components near the interface. With that, there is no dif-ficulty to compute the pressure based on the equation of state, p = " ρE − PNd j=1(ρuj)2 2ρ − γ − bρ γ − 1 B # 1 − bρ γ − 1 . Up to this point, our discussion stresses only on an approach that is capable of keeping the pressure in equi-librium for a model interface-only problem. Since in prac-tice we are interested in shock wave problems as well, we should take the equations, i.e., (12), (13), and (14), in a form such that γ, b, and B remain unchanged across both shocks and rarefaction waves. In this regard, it is known that with γ governed by (13), there is no problem to do so (cf. [1,29]). For b and B, however, due to the appearance of the density term in (12) and (14), it turns out that, in a time when such a situation occurs, for consistent with the mass conservation law of the fluid mixture they should be modified by ∂ ∂t γ − bρ γ − 1 B + Nd X j=1 ∂ ∂xj γ − bρ γ − 1 Buj = γB γ − 1 Nd X j=1 ∂uj ∂xj , (15) and ∂ ∂t bρ γ − 1 + Nd X j=1 ∂ ∂xj bρ γ − 1uj = 0, (16) respectively, so that the mass-conserving property of the solution in the single component region can be acquired also (cf. [30]). We note that, for convenience, we call the set of equations: (13), (15), and (16), a γ-based effective equations for the mixture of the material quantities of the modified Noble-Abel equation of state to be distinct from the other one presented below.
2.2.2Y -based effective equations
Before proceeding further, we comment that to find the initial fluid mixtures, 1/(γ − 1), (γ − bρ)B/(γ − 1), and bρ/(γ −1), that is necessary when we initialize the data for multicomponent flow computations, we use the equation of state (8), where written as a function of the volume fraction Yι, for ι = 1, 2, and Y1+ Y2= 1, it reads
1 − bρ γ − 1p + γ − bρ γ − 1B = ρe = 2 X ι=1 Yιριeι = 2 X ι=1 Yι 1 − bιρι γι− 1 pι+ γι− bιρι γι− 1 B ι .
Here the subscript “ι” denotes the state variable of fluid component ι. By taking a similar approach as employed in Section 2.2.1 for the derivation of the γ-based effec-tive equation it comes out readily a splitting of the above expression into the relations:
1 γ − 1 = 2 X ι=1 Yι γι− 1 , γ − bρ γ − 1B = 2 X ι=1 Yι γι− bιρι γι− 1 Bι , bρ γ − 1 = 2 X ι=1 Yι bιρι γι− 1 , (17)
where in the process of splitting the terms the pressure p is chosen to fulfill the condition
1 − bρ γ − 1 p = 2 X ι=1 Yι 1 − bιρι γι− 1 pι. (18)
Notice that when each of the partial pressures is in an equilibrium state within a grid cell, in conjunction with the first and third parts of (17), the pressure p obtained from (18) would remain in the same equilibrium as well, i.e., p = pι , for ι = 1, 2.
Now if the above volume-fraction notion of the states 1/(γ − 1), (γ − bρ)B/(γ − 1), and bρ/(γ − 1) are being employed in the γ-based effective equations, we are able to rewrite them straightforwardly into a componentwise form as ∂ ∂t Yι γι− 1 + Nd X j=1 uj ∂ ∂xj Yι γι− 1 = 0, (19) ∂ ∂t Yι γι− bιρι γι− 1 Bι + Nd X j=1 ∂ ∂xj Yι γι− bιρι γι− 1 Bι uj = Yι γιBι γι− 1 Nd X j=1 ∂uj ∂xj , (20) ∂ ∂t Yι bιρι γι− 1 + Nd X j=1 ∂ ∂xj Yι bιρι γι− 1 uj = 0, (21)
for ι = 1, 2. Then based on the fact that all the material quantities γι, bι, and Bι will be kept as a constant in each
phase of the domain at all time, from (19), it is easy to find the transport equation for the volume fraction Yι as
∂Yι ∂t + Nd X j=1 uj ∂Yι ∂xj = 0, (22)
whereas from (21), we find the mass conservation law for the fluid component ι as
∂ ∂t(ριYι) + Nd X j=1 ∂ ∂xj (ριYιuj) = 0. (23)
Furthermore, in the case of (20), after some simple alge-braic manipulations, it can be expressed by the form
γιBι γι− 1 ∂Yι ∂t + Nd X j=1 uj ∂Yι ∂xj − bιBι γι− 1 ∂ ∂t(ριYι) + Nd X j=1 ∂ ∂xj (ριYιuj) = 0,
yielding naturally the same set of governing equations as before, i.e., (22) and (23). It is apparent that, if the so-lutions of Yι and ριYι are known from the equations, we
may therefore compute 1/(γ − 1), (γ − bρ)B/(γ − 1), and bρ/(γ−1) directly according to (17). Thus, instead of using the γ-based effective equations, it is a viable alternate to use the Y -based equations: (22) and (23), for the motion of the mixture of the material quantities of the problem.
2.2.3 Complete model system
It is important to mention that to make up a complete model system that is capable of dealing with all the fluid component cases, Y = 1 (barotropic phase), Y = 0 (barotropic phase), or Y ∈ (0, 1) ((barotropic and non-barotropic coexistence phase), we have to know the ap-proximate 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 (22) for each of Y1 and Y2 can be combined,
without affecting anything, to a single one for Y as ∂Y ∂t + Nd X j=1 uj ∂Y ∂xj = 0, (24)
leading to the evolution equation we use in the algorithm for that matter. Note that, in devising a fluid-mixture type algorithm for multicomponent problems, one com-mon practice is to consider the mixture of the total den-sity, ρ = ρ1Y1+ ρ2Y2, as one of the basic variables in the
proposed model system (cf. [29,30,32]). When this is the case, it should be more sensible to include only one of the equations in (23) together with the last two relations
in (17) than the use of both (15) and (16) for the variables (γ − bρ)B/(γ − 1) and bρ/(γ − 1).
Putting all the things together, with the equation of state (9), the model equations with which we propose to solve the present two-fluid flow problems in more than one space dimension take the form
∂ρ ∂t + Nd X j=1 ∂ ∂xj (ρuj) = 0 ∂ ∂t(ρui) + Nd X j=1 ∂ ∂xj (ρuiuj+ pδij) = 0 ∂ ∂t(ρE) + Nd X j=1 ∂ ∂xj (ρEuj+ puj) = 0 if Y 6= 1 ∂ ∂t(ρ1Y ) + Nd X j=1 ∂ ∂xj (ρ1Y uj) = 0 ∂Y ∂t + Nd X j=1 uj ∂Y ∂xj = 0, (25)
for i = 1, 2, · · · , Nd. Note that when Y 6= 1 we have a
sys-tem of Nd+ 4 equations for the motion of both the pure
non-barotropic phase and the fluid mixture. Clearly, the first Nd+2 of them are simply the basic conservation of the
mass, momenta (Ndof them), and total energy, while the
last two equations are (23) and (24) that are used to find the variables 1/(γ − 1), (γ − bρ)B/(γ − 1), and bρ/(γ − 1) from the relations in (17), and also the approximate lo-cation of the interface (the volume fraction function Y does it). On the other instance, when 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 vari-ables of interest, including the pressure from the equation of state. The initialization of the state variables in (25) for fluid-mixture cells can be made in a standard way as described in [29,30] for numerical simulations.
It should be remarked that for an easy extension of the proposed model to more complex equations of state we have taken (25) in a form that is analogous to the five-equation model advocated by Allaire et al. [2] for com-pressible multicomponent flow problems with general non-barotropic equation of state. For any given Nd, if the state
variables of the flow are all in the region of the thermody-namic stability (this is the case we are interested in here), it is not difficult to show that (25) is a hyperbolic system in the sense that any linear combination of the matrices Aj, j = 1, 2, · · · , Nd, appearing in the quasi-linear form of
the equations ∂q ∂t + Nd X j=1 Aj(q) ∂q ∂xj = 0 (26)
has real eigenvalues and a complete set of eigenvectors. As an example, we consider the three-dimensional case Nd = 3, and then have the state vector q in (26) defined
by
q = (ρ, ρu1, ρu2, ρu3, ρE, ρ1Y, Y )T,
and the matrices Aj, for j = 1, 2, 3, defined by
A1= 0 1 0 0 0 0 0 K − u2 1 2u1+ η2 η3 η4 Γ ϕ χ −u1u2 u2 u1 0 0 0 0 −u1u3 u3 0 u1 0 0 0 u1(K − H) H + u1η2u1η3u1η4u1(Γ + 1) u1ϕ u1χ −u1ρ1Y /ρ ρ1Y /ρ 0 0 0 u1 0 0 0 0 0 0 0 u1 , A2= 0 0 1 0 0 0 0 −u1u2 u2 u1 0 0 0 0 K − u2 2 η2 2u2+ η3 η4 Γ ϕ χ −u2u3 0 u3 u2 0 0 0 u2(K − H) u2η2H + u2η3u2η4u2(Γ + 1) u2ϕ u2χ −u2ρ1Y /ρ 0 ρ1Y /ρ 0 0 u2 0 0 0 0 0 0 0 u2 , A3= 0 0 0 1 0 0 0 −u1u3 u3 0 u1 0 0 0 −u2u3 0 u3 u2 0 0 0 K − u2 3 η2 η3 2u3+ η4 Γ ϕ χ u3(K − H) u3η2u3η3H + u3η4u3(Γ + 1) u3ϕ u3χ −u3ρ1Y /ρ 0 0 ρ1Y /ρ 0 u3 0 0 0 0 0 0 0 u3 .
With that, the eigenvalues and the corresponding eigen-vectors of the matrices are:
for matrix A1, ΛA1= diag λ(1)1 , λ (1) 2 , · · · , λ (1) 7 = diag(u1− c, u1, u1+ c, u1, · · · , u1), RA1= r(1)1 , r (1) 2 , · · · , r (1) 7 = 1 1 1 0 0 0 0 u1− c u1 u1+ c 0 0 0 0 u2 u2 u2 1 0 0 0 u3 u3 u3 0 1 0 0 H − u1c K/Γ H + u1c u2u3−ϕ/Γ −χ/Γ ρ1Y /ρ 0 ρ1Y /ρ 0 0 1 0 0 0 0 0 0 0 1 , for matrix A2, ΛA2= diag λ(2)1 , λ (2) 2 , · · · , λ (2) 7 = diag(u2− c, u2, u2+ c, u2, · · · , u2), RA2= r(2)1 , r(2)2 , · · · , r(2)7 = 1 1 1 0 0 0 0 u1 u1 u1 1 0 0 0 u2− c u2 u2+ c 0 0 0 0 u3 u3 u3 0 1 0 0 H − u2c K/Γ H + u2c u1u3−ϕ/Γ −χ/Γ ρ1Y /ρ 0 ρ1Y /ρ 0 0 1 0 0 0 0 0 0 0 1 ,
and for matrix A3,
ΛA3 = diag λ(3)1 , λ (3) 2 , · · · , λ (3) 7 = diag(u3− c, u3, u3+ c, u3, · · · , u3), RA3 = r(3)1 , r(3)2 , · · · , r(3)7 = 1 1 1 0 0 0 0 u1 u1 u1 1 0 0 0 u2 u2 u2 0 1 0 0 u3− c u3 u3+ c 0 0 0 0 H − u3c K/Γ H + u3c u1u2−ϕ/Γ −χ/Γ ρ1Y /ρ 0 ρ1Y /ρ 0 0 1 0 0 0 0 0 0 0 1 ; Ajrk(j) = λ(j)k rk(j), j = 1, 2, 3, and k = 1, 2, · · · , 7. Here
c = p[γ/(1 − bρ)][(p + B)/ρ] is the speed of sound of the fluid, and the other notations appeared in the above formulae are set by Γ = (γ − 1)/(1 − bρ), K = (ΓP3j=1u2
j/2) + ζ, H = E + p/ρ, ηi = −uiΓ ,
for i = 1, 2, 3, ϕ = Γ ∆[bB/(γ − 1)] − p∆[b/(γ − 1)]/Γ , and χ = −Γ ∆[γB/(γ − 1)] + p∆[1/(γ − 1)]/Γ , where ζ = −b2/(γ2− 1)[(Γ B2) + (p/Γ )] and ∆[κ] = κ1− κ2.
For the ease of the latter reference, it is customary to write (25) into a more compact expression by
∂q ∂t + Nd X j=1 fj ∂ ∂xj , q = 0, (27)
where fj is taken as the vector-value function of the form
fj= ∂ ∂xj (ρuj) , ∂ ∂xj (ρu1uj+ pδ1j) , · · · , ∂ ∂xj (ρuNduj +pδNdj) , ∂ ∂xj (ρEuj+ puj) , ∂ ∂xj (ρ1Y uj) , uj ∂Y ∂xj T , for j = 1, 2, · · · , Nd. It is easy to see that the functions fj
defined above reduce to the standard flux functions for a single component flow.
3 Numerical approximations
To find approximate solutions of our two-fluid model pre-sented in Section 2, we use a standard high-resolution wave propagation method developed by LeVeque [16,17] for general hyperbolic systems of partial differential equa-tions. This method is a different form of the fluctuation-and-signal scheme of Roe [26,27] in that we solve one-dimensional Riemann problems at each cell interface, and use the resulting waves (i.e., discontinuities moving at con-stant speeds) to update the solutions in neighboring grid cells. To achieve second-order accurate on smooth solu-tions, and sharp and monotone profiles on discontinuous solutions, we introduce slopes and limiters to the method as in many other high-resolution schemes for conservation laws [10,18,36].
3.1 HLL-type approximate Riemann solvers
As a preliminary, we begin by reviewing the construc-tion of the approximate Riemann problem soluconstruc-tions in one space dimension which is one of the major elements in our numerical algorithm. If we consider the case Nd = 3 as
an example, to simplified the notation, rather than using xj for the spatial variable and uj for the particle
veloc-ity, for j = 1, 2, 3, we take the often-employed symbols (x, y, z) and (u, v, w) for that matter instead. Then the problem to be solved is to find the solution of (25) in the direction normal to one of the xy, yz, and xz planes, with piecewise constant data qL and qRto the left and right of
the cell interface. It should be noted that in the applica-tions 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 con-tact discontinuity; there is no vacuum region occurring in the solution. Without loss of generality, in the following we only look at the Riemann problem in the direction normal to the yz plane.
There are a wide variety of approaches proposed in the literature for determining an approximate solution of the Riemann problem, see [38], for example. Here we are inter-ested in a simple variant of the Riemann solver based on the work of Harten, Lax, and van Leer [12] for hyperbolic systems of conservation laws,
∂q ∂t +
∂
∂xf (q) = 0, (28) where q ∈ lRn 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 λLand λR
to the left and right, λL6= λR, separating three constant
states in the space-time domain. If we assume further that λL and λR are known a priori by some simple estimates
based on the local information of the wave speeds, then by using the integral form of the conservation laws over a sufficiently large control volume [−M, M] × [0, T ], for some M and T ∈ lR, it is easy to find the constant state in the middle region, denoted by qm, as the average of the
exact solution over the interval [T λL, T λR] at time T ,
qm= 1 (λR− λL)T Z λRT λLT q(x, T ) dx = λRqR− λLqL− f(qR) + f (qL) λR− λL , (29)
where f (qι) is the flux evaluated at the state qι, for ι = L
or R.
In spite of the fact that the above 2-wave HLL solver has been applied quite successfully to many numerical methods for solving problems governed by (28), it is known however that the numerical result obtained by using this solver is too diffusive for contact discontinuities. In ad-dition to that, because of the lack of information on the
structure of the material interfaces, it is not feasible at all to the more general multicomponent problems like the one considered in this article. Nevertheless, as recommended by Toro et al. [5,39], this 2-wave Riemann solver can be improved quite easily by introducing an additional middle wave of speed um in the solution structure for modeling
the speed of contact discontinuity, yielding a 3-wave HLL (or called HLLC) solver. Note that if we have had the first two components of qm computed from (29), we may
sim-ply set um = q (2)
m /q(1)m, where qm(i) is the ith component
of the vector qm, see [38] for the various other possible
choices.
With that, our goal next is to find the constant states qmL and qmR in the regions mL and mR to the left and
right of the middle wave, respectively. We do this by mak-ing use of 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 taken effect to the current two-fluid model (25) in the prescribed direction, it is not difficult to show that the result is
q(1)mι= f (1) mι λι− um , qmι(2)= umqmι(1), q (3) mι = vιq(1)mι, q(4) mι= wιq(1)mι, q(5)mι= fmι(5)+ um umfmι(1)− fmι(2) λι− um , qmι(6)= fmι(6) λι− um , qmι(7)= q(7)ι , (30)
where fmι(i) = λιqι(i)−f(i)(qι) represents the ith component
of the vector fmι, for ι = L or R. It is interesting to remark
that qmL(i) and q(i)mRsatisfy the basic consistency condition of the integral form of the conservation laws,
um− λL λR− λL q(i)mL+ λR− um λR− λL qmR(i) = q(i) m,
for i = 1, 2, · · · , 6. By taking λL= min (uL− cL, uR− cR)
and λR = max (uL+ cL, uR+ cR), for instance, we may
then set the speed of the three moving discontinuities by λ1 = λL, λ2 = um, λ3 = λR, and the jumps across each
of them by
W1= qmL− qL, W2= qmR− qmL, W3= qR− qmR,
using the result (30). As usual, wave propagation methods (to be described below) are based on using these propa-gating discontinuity to update the cell averages in the cells neighboring each interface.
3.2 Wave propagation methods 3.2.1 One-dimensional case
To review the basic idea of our underlying integration scheme, it is instructive to look at the simplest Nd = 1
case of our model system, ∂q ∂t + f1 ∂ ∂x, q = 0,
with q and f1 defined as in (25) and (27). For simplicity,
we describe the method on a uniform grid with fixed mesh spacing ∆x in the x-direction, and use a standard finite-volume formulation in which the discretized solution Qn j
approximates the cell average over the grid cell [xj, xj+1]
at time tn, Qnj ≈ 1 ∆x Z xj+1 xj q(x, tn) dx.
The time step from the current time tn to the next tn+1
is denoted by ∆t.
In this numerical setup, a first order accurate version of the method in wave-propagation form is a Godunov-type scheme that can be written as
Qn+1j = Qn j − ∆t ∆x mw X m=1 λ− mWm n j+1+ λ + mWm n j. (31)
Here λm and Wm are solutions of the mth wave family,
for m = 1, 2, · · · , mw, obtained from solving the Riemann
problems at cell interfaces xj and xj+1 with a properly
chosen solver, and λ−
= min (λ, 0), λ+ = max (λ, 0). It
is known that method (31) belongs to a class of upwind schemes, and is stable when the typical CFL (Courant-Friedrichs-Lewy) condition is satisfied (cf. [10,18]). More-over, it is not difficult to show that the method is quasi-conservative in the sense that when applying the method to (25) not only the conservation laws but also the trans-port equations are approximated in a consistent manner by the method, see [30] for the details.
To improve the accuracy of (31) to high resolution, it is a standard approach by first introducing correction waves in a piecewise-linear form with zero mean value, and then propagating each wave over the time step ∆t to update the cell averages it overlaps. Without providing the details in this place (cf. [18] for example), with the corrections, (31) is modified by Qn+1j := Qn+1j −1 2 mw X m=1 |µm| (1 − |µm|) fWm n j+1− |µm| (1 − |µm|) fWm n j, (32)
where µm = λm∆t/∆x, and fWm is a limited value of
Wm obtained by comparing Wm with the corresponding
Wmfrom the neighboring Riemann problem to the left (if
λm> 0) or to the right (if λm< 0).
Note that when utilizing the wave-form representation of the Riemann problem solution, it is a common practice to perform the limiting procedure over each component of the wave via a limiter function Φ (e.g., by using the minmod function Φ(θ) = max(0, min(1, θ)) or some others
as discussed in [36]), and set f
Wmj(i) = Φ
θ(i)mj Wmj(i) with θ (i) mj= WmJ(i) Wmj(i) , J = j − 1 if λmj≥ 0 j + 1 if λmj< 0, (33)
where Wmj(i) is the ith component of Wmj. While the above
approach works in a satisfactory manner for many sin-gle component problems, it was shown in [33] that, for general multicomponent problems, to ensure a consistent approximation of the total energy near the interfaces the third limited component of the 2-wave, fW2j(3), should be
modified so that the desired pressure equilibrium can be maintained there.
Having this in mind, following the previous work done in [33], we make a new definition of fW2j(3) as
f W2j(3):= Φ ¯W2J ¯ W2j ¯ W2j− ∆ bB γ − 1 f W2j(4)+ ∆ γB γ − 1 f W2j(5), where ¯W2ι = W2ι(3) + αW (4) 2ι − βW (5) 2ι , for ι = j and J
(i.e., j − 1 or j + 1 depending on the propagating direc-tion of λ2j); α = ∆[bB/(γ − 1)] and β = ∆[γB/(γ − 1)];
∆[κ] = κ1− κ2. With this revision of the limited 2-wave
on the total energy, it is not difficult to show that for a model interface-only problem (see Section 2.2) we again have the required pressure equilibrium without introduc-ing any spurious oscillations in the approximate solution. Besides that we would typically obtain a better resolution of the result as compared to the first order result, see [33] for the details.
3.2.2 Multidimensional case
To extend the one-dimensional wave propagation method to more space dimensions, here we take a simple dimensional-splitting approach in which a multidimen-sional problem is split into a sequence of one-dimenmultidimen-sional problems. Consider the three-dimensional case Nd = 3,
for example. The hybrid barotropic and non-barotropic two-fluid flow problem modeled by (27), i.e.,
∂q ∂t + f1 ∂ ∂x, q + f2 ∂ ∂y, q + f3 ∂ ∂z, q = 0, can be split into
x-sweeps: ∂q ∂t + f1 ∂ ∂x, q = 0, (34a) y-sweeps: ∂q ∂t + f2 ∂ ∂y, q = 0, (34b) z-sweeps: ∂q ∂t + f3 ∂ ∂z, q = 0. (34c) Assuming a uniform Cartesian grid with fixed mesh spac-ing ∆x, ∆y, and ∆z in the x-, y-, and z-direction, re-spectively. In a finite-volume formulation of the solution,
the value Qn
ijk would approximate the cell average of the
solution over the (i, j, k)th grid cell at time tn,
Qn ijk≈ 1 ∆x∆y∆z Z Z Z Cijk q(x, y, z, tn) dx dy dz,
where Cijk = [xi, xi+1] × [yj, yj+1] × [zk, zk+1] denotes
the cubical region occupied by the grid cell (i, j, k). Then a dimensional-splitting (or called Godunov-splitting) ver-sion of the first-order wave propagation method in three dimensions can be written as
Q∗ ijk= Qnijk− ∆t ∆x mw X m=1 λ(1)− m Wm(1) n i+1,jk+ λ(1)+ m Wm(1) n ijk, (35a) Q∗∗ijk= Q ∗ ijk− ∆t ∆y mw X m=1 λ(2)−m Wm(2) ∗ i,j+1,k+ λ(2)+m Wm(2) ∗ ijk, (35b) Qn+1ijk = Q∗∗ ijk− ∆t ∆z mw X m=1 λ(3)−m Wm(3) ∗∗ ij,k+1+ λ(3)+m Wm(3) ∗∗ ijk. (35c)
Note that in the x-sweeps we start with cell average Qn
ijk at time tn and solve (34a) along each row of cells
Cijk with j and k fixed, updating Qnijk to Q ∗
ijk by the
use of (35a), where λ(1)m,ijk and Wm,ijk(1) are solutions of
the mth wave family obtained from solving the one-dimensional Riemann problems in the direction normal to the cell interface between Cijk and Ci+1,jk with Qnijk
and Qn
i+1,jk as initial data. Consequently in the y-sweeps
we can use the Q∗
ijk values as data for solving (34b)
along each column of cells Cijk with i and k fixed, which
gives us Q∗∗
ijk from (35b) . Finally, in the z-sweeps we
use the Q∗∗
ijk values as data for solving (34c) along the
other column of cells Cijk with i and j fixed, yielding the
solution of the next time step Qn+1ijk from (35c). Clearly, in each one-dimensional sweep, one may apply the same high-resolution approach as before to improve numerical accuracy of this splitting method.
It is a known fact that, except for some simple prob-lems, there will be generally splitting error of the method just described. But from numerical experiences it turns out that the splitting error is often no worse than the er-rors introduced by the numerical methods in each sweep, and hence as a practical matter it is typically not nec-essary to use a more accurate splitting approach such as the Strang splitting [35] to reduce the splitting error, see [18] for some discussion of why one might not want to use a higher order splitting method. In addition, we have also observed good results for many practical prob-lems obtained using the present method as compared to the fully discrete wave propagation method (cf. [15–17]). For these reasons, we will use the one-dimensional high-resolution wave propagation method together with the
Godunov splitting for all the multidimensional tests done in the next section.
4 Numerical results
We now show some numerical results to validate our two-fluid algorithm described in Section 3. In all the problems considered below, we carry out the test using the HLLC Riemann solver, the ”minmod” limiter in a high-resolution method, and the Courant number 0.9.
Example 4.1. Our first example concerns a model shock-contact interaction problem in one dimension that verifies convergence of our numerical solutions to the weak ones in a multicomponent case. Here the initial condition we use is composed of a stationary air-water interface at x = 0.6m within a one meter length domain, and a right-ward going Mach 1.94 shock wave at x = 0.5m traveling from the left to right. The fluid on the right of the interface is an air with the state variables
(ρ, p, Y )R=1.2 kg/m3, 105Pa, 0,
and the material quantities γ = 1.4, b = 2 × 10−3m3/kg,
while the fluid on the left of the interface, (i.e., on the middle and the preshock state), is a water with the state variables
(ρ, p, Y )M =103kg/m3, 105Pa, 1. and the material quantities γ = 7, B = 3 × 108Pa, ρ
0 =
103kg/m3, p
0 = 105Pa. The state behind the shock is
taken to be
(ρ, u, p, Y )L=1337.61 kg/m3, 710.525 m/s, 2 × 109Pa, 1;
see the dashed line shown in Fig. 1 for illustration. We note that in this instance since the interface is accelerated by a shock wave coming from the heavy-fluid to the light-fluid region, it is known that the resulting wave pattern after the interaction would consist of a transmitted shock wave, an interface, and a reflected rarefaction wave.
A snap shot of the computed solutions for ρ1Y , ρ−ρ1Y ,
u, and p are shown in Fig. 1 at time t = 200µs, where we solve the problem using the high-resolution method with 400 mesh points. Observing the displayed profiles, we see clearly the good behavior of the computed air-water interface and also the shock and rarefaction waves as in comparison with the exact solution.
Example 4.2.We are next concerned with a radially symmetric problem where the computed solutions in two space dimensions can be compared to the one-dimensional results for numerical validation. In this test, we use the following two-phase (gas-liquid) flow data for experiments in which, in the gas phase, the state variables are
0.2 0.4 0.6 0.8 0 0.2 0.4 0.6 0.8 1 1.2 ρ 1 Y (Kg/m 3 ) 0.2 0.4 0.6 0.8 0 1 2 3 4 5 6 ρ − ρ 1 Y (kg/m 3 ) 0.2 0.4 0.6 0.8 0 0.2 0.4 0.6 0.8 1 1.2 x (m) u (km/s) 0.2 0.4 0.6 0.8 0 0.5 1 1.5 2 x (m) p (GPa) air water
Fig. 1. High-resolution results for a model shock-contact in-teraction problem at time t = 200µs. The solid line is the exact solution and the points show the computed solution with 400 mesh points. The dashed line in each subplot is the initial con-dition at time t = 0.
and
(ρ, p, Y ) =1.2 kg/m3, 105Pa, 1
if r < r1and r > r2, respectively, while in the liquid phase
they are (ρ, p, Y ) =103kg/m3, 105Pa, 0 if r1< r ≤ r2, where r = p x2+ y2, r 1= 0.2m, and r2=
0.7m. Here the material quantities of the fluid is taken the same as before, with the exception that a smaller value of the constant for the excluded volume, b = 10−4m3/kg, is
used for the gas.
We note that for this problem all the fluid component is in a resting state initially with zero total velocity, but due to the pressure difference between the fluids at r = r1,
breaking of the inner circular membrane occurs instanta-neously, yielding an outward-going shock wave in liquid, an inward-going rarefaction wave in gas, and a contact discontinuity lying in between that separates the gas and liquid. As times go along, the inward-going wave would be reflected from the geometric center that generates a new outward-going wave and induces the subsequent interac-tion of waves. At a somewhat later time, the outward-going shock wave would be collided with the outer gas-liquid interface at r = r2 that results in a wave pattern
consisting of a transmitted shock wave, an interface, and a reflected rarefaction wave. Because of the symmetry of the solution, for simplicity, we only take a quarter of the unit square, and make use of the line of symmetry bound-ary conditions to the bottom and the left sides during the computations.
Figures 2 and 3 show numerical results for the total density, radial velocity (defined as ¯u = √u2+ v2), and
pressure, at three stopping times, t = 120, 240, and 450µs, where the test has been carried out by using a 200 × 200
grid with the high-resolution method. Clearly, from the contour plots shown in Fig. 2, we observe good resolution of the solution structure (i.e., both the shock and interface remain circular and appear to be very well located) after the breaking of the membrane and also the interaction of the shock and the outer interface.
The scatter plots shown in Fig. 3 provide the validation of our two-dimensional results as in comparison with the “true” solution obtained from solving the one-dimensional model with appropriate source terms for the radial sym-metry, using the high-resolution method with 1000 mesh points in a unit length domain. That is, for the equation, we have a modified-version of the model (27) in one di-mension as ∂q ∂t + f ∂ ∂r, q = ψ(q) (36) with f a vector-value function defined by
f = ∂ ∂r(ρu) , ∂ ∂r ρu 2+ p, ∂ ∂r(ρEu + pu) , ∂ ∂r(ρ1Y u) , u ∂Y ∂r T ,
and ψ the source term derived directly from the geomet-ric simplification of a multidimensional flow to a one-dimensional one,
ψ = −αr ρu, ρu2, ρEu + pu, ρ1Y u, 0
T
. Note that in the case of a 2D radially or 3D spherically symmetric flow, we use the quantity α = 1 or 2, re-spectively; u now denotes the particle velocity in the r-(radial) direction. We use a Strang-type time splitting pro-cedure [35] to deal with the geometric sources of (36) in a high-resolution manner during the run. From the figure, it is clear that our results agree quite well with the ”true” solutions at all the selected times, and also free of wrong fluctuations in the pressure near the inner and outer inter-faces before and after the interactions of rarefaction and shock waves.
Example 4.3.To show how our algorithm performs on shock waves in a more general two-dimensional geometry, we are interested in the simulation of a shock-induced col-lapse of a cylindrical air cavity in water studied by Bourne and Field [6] experimentally, and Ball et al. [4] numer-ically. We setup the problem by introducing a station-ary air cavity of radius r0 = 3mm located at (x0, y0) =
(6, 0)mm and also a planarly rightward-going shock wave in water at x = 1mm approaching toward the cavity on the left. Similar to the initial condition employed in Exam-ple 4.1, inside the air bubble, we use the data on the state R, while outside the air bubble (the fluid is water), we use the data on the states M and L, when they are in the pre-shock and post-pre-shock regions, respectively. Here, in car-rying out the computations below, we take a shock tube of size [0, 15] × [0, 6](mm)2, and use the solid wall
bound-ary condition on the bottom side, and the non-reflecting boundary condition on the remaining sides.
Density
t=120µ s
Radial velocity Pressure
t=240µ s
t=450µ s
Fig. 2. High-resolution results for a two-phase radially sym-metric problem. Contour plots for the total density, radial ve-locity, and pressure are shown at three different times t = 120, 240, and 450µs. The dashed line shown in the graph is the approximate location of the gas-liquid interface.
0.2 0.4 0.6 0.8 0.2 0.4 0.6 0.8 1 1.2 ρ (Kg/m 3 ) t= 120 µ s 0.2 0.4 0.6 0.8 t= 240 µ s 0.2 0.4 0.6 0.8 t= 450 µ s
gas liquid gas
0.2 0.4 0.6 0.8 0 0.1 0.2 0.3 0.4 u (km/s) 0.2 0.4 0.6 0.8 0.2 0.4 0.6 0.8 0.2 0.4 0.6 0.8 0 0.2 0.4 0.6 0.8 1 p (GPa) r (m) 0.2 0.4 0.6 0.8 r (m) 0.2 0.4 0.6 0.8 r (m)
Fig. 3.Scatter plots of the results for the run shown in Fig. 2. The solid line is the “true” solution obtained from solving the one-dimensional model with appropriate source terms for the radial symmetry using the high-resolution method and 1000 mesh points. The dotted points are the two-dimensional result. The dashed line is the approximate location of the gas-liquid interface.
Figures 4 shows the Schlieren image of the results for the density and pressure at ten different times t = i × 0.4µs, for i = 1, 2, · · · , 10 (measured relative to a time when the incident shock wave first hits the upstream wall of the air cavity), using the same high-resolution method as in the previous example and a 1500×600 grid. From the figure, it is easy to see that after the passage of the shock to the cavity, the upstream wall begins to spall across the cavity, yielding an refracted air shock traveling within it at the speed Vr ≈ 2.231km/s until its first reflection
on the far cavity wall at the time t ≈ 2.647µs (see the comments given below on how to obtain these numerical values). Noticing that this upstream cavity wall would in-volute eventually to form a jet which subsequently crosses the cavity and sends an intense blast wave out into the sur-rounding liquid upon impacting with the far cavity wall. Following from [6], the time from the initial impingement of the shock on the upstream wall to the impact of the jet on the downstream wall is referred to as the collapsed time, and in the current case it is about 2.787µs. At a later stage of the collapse, highly compressed air would be trapped in a lobe-like region (see the plots at times t = 3.6 and 4µs), where high pressure and also high tem-perature are concentrated on. It should be mentioned that this lobe-like region is the place where light emission was observed in an experimental work of Bourne and Field [6]. As the time goes along, complex multiple reflection of waves would be continued to happen. The cross section of the results for the same run along the bottom boundary is drawn in Fig. 5 giving some quantitative information about the density and pressure at the selected times. As far as the global wave structure of the solution is con-cerned the results presented here are reasonable one as compared to the one appearing in Ball et al. [4], where a Free-Lagrange method was employed to solve a similar problem but with a γ-law gas for the air and an incident shock pressure 1.9GPa.
To get a further assessment of the key characteristics of the solutions in the problem, in Fig. 6, we report a diagnosis of the space-time locations of the incident shock wave, the upstream cavity wall, the downstream cavity wall, the first refracted shock wave, and the jet, where we again use a 1500×600 grid and the same numerical method as before. With that, it is a common practice (cf. [21,25]) to perform a linear least-squares fit of these trajectories separately, and take the slope of the respective line as one measure of the speed of Vs, Vu, Vd, Vr, and Vj, in the order
mentioned. Here, at some sample times, we obtain the positions of the incident and refracted shocks by checking the jumps in the pressure to see whether there are greater than a prescribed tolerance, say εp = (pL − pM)/10, at
the horizontal sections of 0.01mm below the top boundary, and the axis of symmetry, respectively. As to the positions of the other features, we check the volume fractions to see if there are close to 1/2 at the axis of symmetry, with the exception that, for the upstream wall, we use a section at a height 1.2mm above the axis when time t ≥ 1.5µs.
For this air-cavity collapse problem, it is of great in-terests to see how the time for the cavity collapse (defined
Fig. 4. Two-dimensional results for a planar shock wave in water over a cylindrical air cavity. Schlieren-type image for the density and pressure are shown at ten different times t = i × 0.4µs, for i = 1, 2, · · · , 10, where a 1500 × 600 grid was used in the computation. The dashed line is the initial location of the air-water interface.
0 5 10 15 −1 0 1 2 3 4 0 0.5 1 1.5 2 t (µs) x (mm) ρ (Kg/m 3 ) 0 5 10 15 −1 0 1 2 3 4 0 2 4 6 8 10 t (µs) x (mm) p (GPa)
Fig. 5. Cross-sectional plots of the results for the run shown in Fig. 4 along the bottom boundary.
0 2 4 6 8 10 12 0 0.5 1 1.5 2 2.5 3 x (mm) t ( µ s) V s V u V d V r V j
Fig. 6. Space-time locations of the incident shock wave, the upstream cavity wall, the downstream cavity wall, the first refracted shock wave, and the jet for the run shown in Fig. 4. These trajectories can be used to estimate the speeds such as Vs, Vu, Vd, Vr, and Vj in the order mentioned above.
before) as well as the velocity of the jet are varied as we change the magnitude of the incident shock pressure. Here we carry out one such computation with the set of the pressures pL = i × 0.5GPa, for i = 1, 2, · · · , 8, and show
the result in Table 1, observing a monotonic decreasing and so increasing of the collapsed time and the jet veloc-ity, respectively, as the incident shock pressure increases. Figure 7 is the graph about the jet velocity together with the Hugoniot locus for the water with the same pre-shock condition as before. Notice that the Hugoniot now divides the graph into two parts where in the lower part the jet velocity is less than the shock velocity and so the incident shock crosses the cavity ahead of the jet, while in the up-per part the jet arrives at the downstream cavity wall first and may send out a shock in the material downstream of the cavity ahead of the incident shock. In case the incident shock pressure is in the close neighborhood of 3.5GPa, the speed of the incident shock would be in the proximity to that of the jet. A figure of this kind has also been displayed in the experimental work of Bourne and Field [6], giving some indication of the sensible behavior of the solutions on the computed jet velocity. It should be remarked that Bagabir and Drikakis [3] have done a similar numerical
Table 1. Data summarizing solutions of a shock-induced air cavity problem for the particle velocity, jet velocity, and col-lapsed time, at incident shock pressures pL= i × 0.5GPa, for
i= 1, 2, · · · , 8
Shock pressure Particle velocity Jet velocity Collapsed time (GPa) (km/s) (km/s) (µs) 0.5 0.256 1.309 6.091 1.0 0.435 1.817 4.119 1.5 0.582 2.191 3.290 2.0 0.711 2.572 2.787 2.5 0.826 2.729 2.487 3.0 0.933 3.064 2.258 3.5 1.032 3.350 2.072 4.0 1.125 3.625 1.952 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 1 1.5 2 2.5 3 3.5 4 p (GPa) V j (km/s)
Fig. 7. The variation of the jet velocity, denoted by Vj for
the shock-induced air-cavity problem, with incident shock pres-sures pL= i × 0.5GPa, for i = 1, 2, · · · , 8. The solid line shown
in the graph is the Hugoniot locus for water with same pre-shock condition as in Fig. 4, and the dotted points are the numerical results.
study, but used the Mach number as the basic parameter for a shock wave in air over a cylindrical helium cavity.
Finally, Fig. 8 shows numerical results for a study of the time history of the volume fraction of the air cavity at four different incident shock pressures pL = 1, 2, 3, and
4GPa. As a close look at the data for each of the displayed volume fraction, we may find a time interval during which Y falls from 0.8 to 0.1 at almost a linear rate. Consider the case when PL= 2GPa, for example, the time interval
for that to occur is t ∈ [0.862, 2.933]µs. Our results also show that not long after the collapsed time the volume continuing to decline, but at a somewhat reduced rate which is in consistent with many results shown for shock-driven collapse of cylindrical cavities in fluid dynamical problems (cf. [4,6,11]). A detailed study of this problem for a spherical air cavity case in three dimensions will be reported elsewhere.
5 Conclusions
We have described an Eulerian interface-capturing ap-proach for a hybrid barotropic and non-barotropic two-fluid flow problem in more than one space dimension. The algorithm uses a volume-fraction formulation of the model equations together with an extended equation of state that is devised to ensure a consistent approxima-tion of the energy equaapproxima-tion near the numerical-induced
0 1 2 3 4 5 6 0 0.2 0.4 0.6 0.8 1 t (µs) Y p=1GPa p=2GPa p=3GPa p=4GPa
Fig. 8. The time history of the volume fraction of the air cavity at four different incident shock pressures pL= 1 ,2 ,3,
and 4GPa.
smeared material interfaces, and an easy computation of the pressure from the equation of state. A high-order Go-dunov method based on the wave-propagation viewpoint is employed to solve the proposed model system with the dimensional-splitting technique included in the method for multidimensional problems. Numerical results presented in this paper demonstrate clearly the feasibility of the ap-proach for a reasonable class of two-fluid problems with the thermodynamic property of the barotropic and non-barotropic fluid components characterized by the Tait and Noble-Abel equations of state, respectively.
Acknowledgement. This work was supported in part by the National Science Council of Republic of China Grants NSC-89-2115-M-002-033, 90-2115-M-002-022, and 91-2115-M-002-016.
References
R. Abgrall. How to prevent pressure oscillations in multicom-ponent flow calculations: a quasi conservative approach. J. Comput. Phys., 125:150–160, 1996.
G. Allaire, S. Clerc, and S. Kokh. A five-equation model for the simulation of interface between compressible fluids. J. Comput. Phys., 181:577–616, 2002.
A. Bagabir and D. Drikakis. Mach number effects on shock-bubble interaction. Shock Waves, 11:209–218, 2001. G. J. Ball, B. P. Howell, T. G. Leighton, and M. J. Schofield.
Shock-induced collapse of a cylindrical air cavity in water: a Free-Lagrange simulation. Shock Waves, 10:265–276, 2000. P. Batten, N. Clarke, C. Lambert, and D. Causon. On the choice of wave speeds for the HLLC Riemann solvers. SIAM J. Sci. Comput., 18:1553–1570, 1997.
N. K. Bourne and J. E. Field. Shock-induced collapse of single cavities in liquid. J. Fluid Mech., 244:225–240, 1992. J. P. Dear and J. E. Field. A study of the collapse of arrays of
cavities. J. Fluid Mech., 190:409–425, 1988.
R. P. Fedkiw, T. Aslam, B. Merriman, and S. Osher. A non-oscillatory Eulerian approach to interfaces in multimaterial flows (the ghost fluid method). J. Comput. Phys., 152:457– 492, 1999.
E. Fermi. Thermodynamics. Dover, New York, 1956.
E. Godlewski and P.-A. Raviart. Numerical Approximation of Hyperbolic Systems of Conservation Laws. Springer-Verlag, 1996. Applied Mathematical Science 118.
J.-F. Haas and B. Sturtevant. Interaction of weak shock waves with cylindrical and spherical gas inhomogeneities. J. Fluid Mech., 181:41–76, 1987.
A. Harten, P. D. Lax, and B. van Leer. On upstream differenc-ing and Godunov-type schemes for hyperbolic conservation laws. SIAM Review, 25:35–61, 1983.
S. Karni. Multicomponent flow calculations by a consistent primitive algorithm. J. Comput. Phys., 112:31–43, 1994. B. Koren, M. R. Lewis, E.H. van Brummelen, and B. van
Leer. Riemann-problem and level-set approaches for ho-mentropic two-fluid flow computations. J. Comput. Phys., 185:654–674, 2002.
J. O. Langseth and R. J. LeVeque. A wave propagation method for three-dimensional hyperbolic conservation laws. J. Comput. Phys., 165:126–166, 2000.
R. J. LeVeque. High resolution finite volume methods on arbi-trary grids via wave propagation. J. Comput. Phys., 78:36– 63, 1988.
R. J. LeVeque. Wave propagation algorithms for multi-dimensional hyperbolic systems. J. Comput. Phys., 131:327–353, 1997.
R. J. LeVeque. Finite Volume Methods for Hyperbolic Prob-lems. Cambridge University Press, 2002.
D. R. Lide. Handbook of Chemistry and Physics. CRC Press, 76 edition, 1996.
H. Luo, J. D. Baum, and R. L¨ohner. On the computation of multi-material flows using ALE formulation. J. Comput. Phys., 194:304–328, 2004.
A. Marquina and P. Mulet. A flux-split algorithm applied to conservative models for multicomponent compressible flows. J. Comput. Phys., 185:120–138, 2003.
S. P. Marsh. LASL Shock Hugoniot Data. University of Cali-fornia Press, Berkeley, 1980.
F. D. Murnaghan. Finite Deformation of an Elastic Solid. John Wiley & Sons, Inc., New York, 1951.
H. Nagoya, T. Obara, and K. Takayama. Underwater shock wave propagation and focusing in inhomogeneous media. In R. Brun and L. Z. Dumitrescu, editors, Proceedings of 19th Intl. Symp. on Shock Waves, Marseille, pages 439– 444. Springer-Verlag, Berlin, 1995.
J. J. Quirk and S. Karni. On the dynamics of a shock-bubble interaction. J. Fluid Mech., 318:129–163, 1996.
P. L. Roe. Fluctuations and signals– A framework for nu-merical evolution problems. In K. W. Morton and M. J. Baines, editors, Numerical Methods for Fluid Dynamics, pages 219–257. Academic Press, 1982.
P. L. Roe. Upwind schemes using various formulations of the Euler equations. In F. Angrand, A. Dervieux, J. A. Desideri, and R. Glowinski, editors, Numerical Methods for the Euler Equations of Fluid Dynamics, pages 14–31. SIAM, 1985.
R. Saurel and R. Abgrall. A simple method for compress-ible multifluid flows. SIAM J. Sci. Comput., 21:1115–1145, 1999.
K.-M. Shyue. An efficient shock-capturing algorithm for com-pressible multicomponent problems. J. Comput. Phys., 142:208–242, 1998.
K.-M. Shyue. A fluid-mixture type algorithm for compressible multicomponent flow with van der Waals equation of state. J. Comput. Phys., 156:43–88, 1999.
K.-M. Shyue. A volume-of-fluid type algorithm for compress-ible two-phase flows. In Hyperbolic Problems: Theory, Nu-merics, Applications, pages 895–904. Birkh¨auser-Verlag, 1999. Intl. Series of Numerical Mathematics, Vol. 130. K.-M. Shyue. A fluid-mixture type algorithm for
compress-ible multicomponent flow with Mie-Grueneisen equation of state. J. Comput. Phys., 171:678–707, 2001.
K.-M. Shyue. A fluid-mixture type algorithm for barotropic two-fluid flow problems. J. Comput. Phys., 200:718–748, 2004.
R. W. Smith. AUSM (ALE): a geometrically conservative ar-bitrary Lagrangian-Eulerian flux splitting scheme. J. Com-put. Phys., 150:268–286, 1999.
G. Strang. On the construction and comparison of difference schemes. SIAM J. Numer. Anal., 5:506–517, 1968. P. K. Sweby. High resolution schemes using flux limiters for
hy-perbolic conservation laws. SIAM J. Numer Anal., 21:995– 1011, 1984.
E. F. Toro. A fast Riemann solver with constant covolume ap-plied to the random choice method. Intl. J. Numer. Meth. in Fluids, 9:1145–1164, 1989.
E. F. Toro. Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction, 2nd Edition. Springer-Verlag, 1999.
E. F. Toro, M. Spruce, and W. Speares. Restoration of the contact surface in the HLL-Riemann solver. Shock Waves, 4:25–34, 1994.
E.H. van Brummelen and B. Koren. A pressure-invariant con-servative Godunov-type method for barotropic two-fluid flows. J. Comput. Phys., 185:289–308, 2003.
A. B. Wardlaw, Jr. and H. U. Mair. Spherical solutions of an underwater explosion bubble. Shock and Vibration, 5:89– 102, 1998.
K. Yamada, H. Nagoya, and K. Takayama. Shock wave reflection and refraction over a two-fluid interface. In R. Brun and L. Z. Dumitrescu, editors, Proceedings of 19th Intl. Symp. on Shock Waves, Marseille, pages 299–304. Springer-Verlag, Berlin, 1995.
