Article ID jcph.1999.6349, available online at http://www.idealibrary.com on
A Fluid-Mixture Type Algorithm for Compressible Multicomponent Flow with van der Waals
Equation of State
Keh-Ming Shyue
Department of Mathematics, National Taiwan University, Taipei, Taiwan 106, Republic of China E-mail: shyue@math.ntu.edu.tw
Received September 17, 1998; revised May 19, 1999
In previous work by the author, a simple interface-capturing approach has been developed and validated for compressible multicomponent flows with a stiffened gas equation of state in multiple space dimensions. The algorithm uses a mixture type of the model equations written in a quasi-conservative form to ensure a consistent approximation of the energy equation near the interfaces where two or more fluid components are present in a grid cell. A standard high-resolution wave propagation method is employed to solve the proposed system, giving an efficient implementation of the algorithm. In this paper, the method is extended to a more general two-phase (liquid–gas) flow where the fluid of interests is characterized by a van der Waals- type equation of state. Several numerical results are presented in both one and two space dimensions that show the feasibility of the method with the Roe solver as applied to practical problems without introducing any spurious oscillations in the pressure near the interfaces. This includes a convergence study of a shock wave in liquid over a gas bubble. To deal with a difficult slip line problem where there is a strong shear flow moving along the interface, we implement the method based on the shock-only Riemann solver with an additional update by the scheme to the total kinetic energy. Rather than using solutions from the basic conservation laws for the density and momenta which incurs large errors, the resulting total kinetic energy is used to the computation of the pressure from the equation of state, yielding typically more accurate results than the unmodified method near the slip lines. This is demonstrated by numerical results of some sample two-dimensional Riemann problems. ° 1999 Academic Pressc
Key Words: fluid-mixture algorithm; interface capturing method; multicomponent flow; slip line problems; stiffened gas equation of state; van der Waals gases.
43
0021-9991/99 $30.00 Copyright c° 1999 by Academic Press All rights of reproduction in any form reserved.
1. INTRODUCTION
This paper is concerned with the development of a simple interface-capturing approach for efficient numerical resolution of compressible multicomponent flows with a van der Waals-type fluid (i.e., a fluid with the finite size of a molecule and the nonzero cohesive forces between molecules [21]). We consider a simplified two-phase flow where the fluids of interest consist of two different phases, liquid and gas, separated by immiscible interfaces;
see Fig. 1 for a sample setup in two space dimensions. The algorithm uses a Eulerian formulation of the equations in which, on the gas-phase part of the domain, the fluid is governed by the full set of the compressible Euler equations. In two space dimensions, for instance, it takes the form
∂t
ρ ρu ρv ρE
+ ∂x
ρu ρu2+ p
ρuv ρEu + pu
+ ∂y
ρv ρuv ρv2+ p ρEv + pv
= 0, (1)
whereρ is the density, u and v are the particle velocities in the x- and y-direction, respec- tively, p is the pressure, and E is the specific total energy. We assume that the gas satisfies a van der Waals equation of state,
p(ρ, e) =
µγ − 1 1− bρ
¶
(ρe + aρ2) − aρ2, (2)
so as to deal with the possible real-gas effect (without phase transitions) when both the temperature and pressure are high (cf. [21, 74, 78]). Here e denotes the specific internal energy,γ is the ratio of specific heats (γ > 1), and the quantities a, b are the van der Waals gas constants for molecular cohesive forces and the finite size of molecules, respectively (a≥ 0, 0 ≤ b < 1/ρ, see [44] for numerical values to various gaseous substances). As usual we set E= e + (u2+v2)/2. Note that a van der Waals gas of the form (2) reduces to a Noble–
Abel gas (also called a constant covolume gas) when a= 0 [54, 73] and to a polytropic gas when b= 0 as well. The four components of (1) express the conservation of mass, momenta in the x- and y-direction, and energy, respectively [15].
FIG. 1. A typical example of a two-phase flow setup with interfaces that separate regions of two different fluids, liquid and gas, into two parts. Note that the gas component in each elliptical-like shape of the domain may be different from one another. It is those grid cells that are cut by the interfaces requiring special attention for proper numerical treatments.
On the liquid-phase part of the domain, however, while the motion of a liquid is assumed to be governed by (1), the algorithm uses the stiffened gas equation of state,
p(ρ, e) = (γ − 1)ρe − γ B, (3)
for a fundamental characterization of material properties of the liquids. HereB is a pressure- like constant that, together withγ , can be determined by a fitting procedure from laboratory data (cf. [47]). A typical set of parameter values is for water:γ = 7, B = 3000 atm [15, 25], and for human blood: γ = 5.527, B = 614.6 MPa [52], approximately. It should be mentioned that in addition to the modeling of a liquid, Eq. (3) is often used to describe other type of materials, including many compressible solids of practical importance (cf. [47, 55, 60]).
We want to use a state-of-the-art shock-capturing method on a uniform rectangular grid for the multicomponent flow computations. Clearly when grid cells contain only a single phase of the fluid, there is no problem to solve each phase of the equations separately. But in practice due to the presence and the subsequent dynamic-evolution of the interfaces by the solutions of the governing system, it is inevitable to have two or more fluid components staying within a cell (see Fig. 1). Because of this, the need to give a proper modeling and approximation of these mixture grid cells becomes a principle issue in many of the multicomponent algorithms developed in the literature (cf. [1, 12, 20, 29, 30, 46, 60, 63]).
See [60] in particular for a concise survey of the up-to-date multicomponent methods.
The approach we take here is an extension of the work described by Saurel and Abgrall [60] and by the author [63, 64], in that as opposed to a simplier case with the stiffened gas equation of state (3), a modified van der Waals equation of state (5), see Section 2, is introduced as a basic element to the modeling of the mixing between the stiffened and van der Waals gases. With that, assuming uniform pressure equilibrium and constant particle velocity across the interfaces, from the energy equation, we are able to derive the effective equations for the mixture of the material-dependent quantities near the interfaces. As in the previous work [63], we take these equations to be of the form that do not vary their solutions across the shocks and rarefaction waves as well. Combining the resulting set of effective equations to the Euler equations yields a model system that is written in a quasi- conservative form; see (13) and (36) for the one- and two-dimensional models, respectively.
We use the high-resolution wave propagation method developed by LeVeque [34, 39, 43] to solve the proposed system. Numerical results present in Sections 3.4 and 4.3 show that this is a viable approach in both one and two dimensions as the method is applied with the Roe solver to practical problems without introducing any spurious oscillations in the pressure near the interfaces.
To deal with a difficult slip line problem where there is a strong shear flow moving along the interface, we implement the two-dimensional method based on a shock-only Riemann solver with an additional update by the scheme to the total kinetic energy. Rather than using the solutions from the basic conservation laws for the density and momenta, which incurs large errors, see Section 5.1, the resulting total kinetic energy is employed to compute the pressure from the equation of state, yielding typically more accurate results than the unmodified method near slip lines. This will be discussed further in Section 5.2.
It is true that in many multicomponent flow applications physical effects, such as surface tension and viscosity, play an important role to not only controlling the dynamics of the interfaces, but also influencing the structure of the nearby solutions, and hence need to be
taken in the model for realization (cf. [66]). Typical examples of this kind are the popular water-drop problem in air and the rising air bubble problem in water. In these instances, because the fluids are mostly in a low Mach number regime, compressibility of the fluids is often ignored, and so the problem may be formulated into an incompressible form and solved by a state-of-the-art numerical method for incompressible flows; see [7, 68, 69, 75]
for an example.
Here in contrast to the work just mentioned, we consider a class of problems where the influence of compressibility of the fluids to the solutions is vital, but not the surface tension and viscosity. Examples cover a family of shock wave problems with complicated interface patterns [24, 25, 53, 71, 76, 77] and a hydrodynamic model of sonoluminescence (an acoustic-induced light emission phenomenon) [8, 33, 78]. It is interesting to note that the latter problem in fact motivates the current study of two-phase flows with a van der Waals-type fluid where (2) is used for the modeling of real gases and (3) is employed for a simple approximation of liquids. Our goal here is to establish a basic solution strategy and validate its use via numerical experimentation of some sample problems. Direct simulation of sonoluminescence and other important two-phase flow problems such as shock waves in bubbly liquids or liquid–solid suspensions will be considered in the future; see [65] for a preliminary result of the latter problem.
This paper is organized as follows. In Section 2, we begin by discussing the basic equation of state and the associated thermodynamic stability for the mixing of the stiffened and van der Waals gases within a grid cell. In Section 3, we describe the one-dimensional version of the multicomponent algorithm in more details. This includes the construction of the Riemann problem solution by the shock-only or Roe solver (see Section 3.2) and some numerical results that validate the proposed approach (see Section 3.4). Extension of the basic approach to two space dimensions is explained briefly in Section 4, and some sample two-dimensional results are shown in Section 4.3. A study of how the algorithm works for slip line problems is addressed in Section 5.1, and a correction of the algorithm is made in Section 5.2.
2. EQUATIONS OF STATE
To begin, we introduce a hybrid version of the equation of state that is necessary in the algorithm for modeling the mixing between the stiffened and van der Waals gases within a grid cell. We do this by taking an approach that expresses (2) and (3) in terms of their natural variables: the specific entropy S and the specific volume V= 1/ρ, yielding the formulas
p(V, S) = A(S)(V − b)−γ − aV−2, (4a) and
p(V, S) = A(S)V−γ − B, (4b)
in a respective manner, whereA(S) = R exp[(S − Sref)/CV] takes the same form in both cases (cf. [21, 26]). Here CV represents the specific heat at constant volume (see below),R is the universal gas constant, and Srefis the reference state of the specific entropy. Assuming the fluids under consideration are all in an adiabatic equilibrium with the same entropy, it is feasible then to define an extension of the equation of state that combines (4a) with (4b)
in the following way:
p(V, S) = A(S)(V − b)−γ − (B + aV−2). (4c) With this, it is easy to see that not only the limiting van der Waals gas case (4a) is recovered as B → 0, but so also is the stiffened gas case (4b) as both a → 0 and b → 0. In the nonlimiting case, however, where none of the material-dependent quantities a, b, and B is close to zero, (4c) does give a way to the representation of the cases in between, that is to the mixing of the stiffened and van der Waals gases.
It should be noted that, as far as the computational efficiency of the method is concerned, (4c) is of little use in practice, for it requires some extra work to evaluate the specific entropy from the solutions of the Euler equations. Instead, by employing the first and second laws of thermodynamics, we may rewrite (4c) in terms of the often-used variables in gas dynamics, ρ and e, as
p(ρ, e) =
µγ − 1 1− bρ
¶
(ρe − B + aρ2) − (B + aρ2). (5)
It is clear that (5) is a generalization of the equations of state (2) and (3). As we will see in the latter sections (cf. [64, 65] also for a simpler case), with this so-called modified van der Waals equation of state (5), it is very robust to devise an interface-capturing solver for a family of two-phase flow problems considered here.
As to the computation of the fluid temperature T (which is important to sonolumines- cence, for example), with (5), we may simply use one of the formulas,
p(V, T ) = RT
V− b − B − a V2, e(V, T ) = RT
γ − 1+ BV − a V,
for realization. Note that these two equations are easily derived from (5) using the basic thermodynamic principles; see [21] for more details.
It is important to mention that, in this work, the thermodynamic description of the multi- component flows is limited by the stability requirement that the internal energy defined in (5) be a convex function of its dependent variables V and S. Because of this, the immediate consequence is the exclusion of the important but difficult problems involving the transition of phases here. Analogously in [49, 50] for a Mie–Gr¨uneisen-type equation of state (see below), it can be shown that explicit conditions for the above mentioned thermodynamic stability of a van der Waals gas are
(i) The heat capacity coefficients CV and Cpmust be positive,
CV = (∂Te)|V = R γ − 1> 0,
Cp = (∂Th)|p = R[γ RT − 2aρ(1 − bρ)2] (γ − 1)[RT − 2aρ(1 − bρ)2] > 0, where h= e + (p/ρ) is the specific enthalpy.
(ii) The isentropic bulk modulus KSmust be positive,
KS= ρ(∂ρp)|S= µ γ
1− bρ
¶
(p + B + aρ2) − 2aρ2> 0.
(iii) The product of the thermal expansion coefficientϑ and the Gr¨uneisen coefficient 0 must be nonnegative,
ϑ0 =
"
−1 ρ(∂Tρ)¯¯
¯¯p
#"
1 ρ(∂ep)¯¯
¯¯ρ
#
=
· R(1 − bρ) RT − 2aρ(1 − bρ)2
¸·γ − 1 1− bρ
¸
= R(γ − 1)
RT − 2aρ(1 − bρ)2 ≥ 0.
In summary, for stability, from the above conditions, this amounts to the satisfaction of the inequalities
γ > 1, RT > 2aρ(1 − bρ)2.
Combining this with the positiveness of the basic thermodynamic states,ρ, T , and 0, defines the domain of the phase space of this van der Waals gas model. Note in particular that from item (ii) given above we have KS= ρc2, and so the positive of KSimplies that the speed of sound c belongs to a set of real numbers.
For convenience, we write (5) into a more general Mie–Gr¨uneisen-type equation of state of the form
p(ρ, e) = 0(ρ)[ρe − ρeH(ρ)] − pH(ρ), (6) where0 is the Gr¨uneisen coefficient of the material of interests, and the density-dependent functions eH, pH are the reference Hugoniot states. Typically, the specific form of (6) is rather complicated in general (cf. [27, 28, 48]). Here for fluids described by (5), we have a relatively simple but nontrivial model to work with, i.e.,
0(ρ) = γ − 1
1− bρ, eH(ρ) = B − aρ2
ρ , pH(ρ) = B + aρ2. (7) It is without question that experience gained by studying the current van der Waals gas model will help to the further development of an efficient multicomponent solver for more general materials as described by (6).
3. ONE SPACE DIMENSION
To motivate the basic idea of our method in multiple space dimensions, we begin by considering one-dimensional problems with both the gas and the liquid governed by the Euler Eqs. (1) of the following form in the x-direction,
∂t
ρ ρu ρE
+ ∂x
ρu ρu2+ p ρEu + pu
= 0, (8)
and by the equations of state (2) and (3), respectively. The algorithm uses a popular approach that employs (8) for the motion of the liquid–gas mixtures of the conserved variablesρ, ρu, and ρE in a multicomponent grid cell (see [1, 61, 79] for the use of other governing equations). We compute the pressure based on the equation of state (6), so long as the mixture of the problem-dependent material quantities appearing in (7), i.e.,γ, a, b, and B, are defined and known a priori as well. In the algorithm, by following a general procedure proposed in [63], conditions for those material quantities are found so that not only the pressure is retained in equilibrium for an interface only fluid-mixture cell, but also the mixture of the total mass remains conservative on the domain (of course, the mass of each fluid component may not typically be conserved here). For completeness, we next give a detailed description of that procedure.
3.1. Derivation of Model Equations
As in the previous works [63, 64], our starting point is to consider an interface only problem where both the pressure p and particle velocity u are constants in the domain (i.e., p and u satisfy the standard surface-tension free dynamic and kinematic boundary conditions on the interface [32], respectively), while the other variables such as the density ρ and the material-dependent parameters in the Mie–Gr¨uneisen equation of state (6) are having jumps across some interfaces. In this setup, we first write (8) in the following nonconservative form,
∂tρ + u∂xρ + ρ∂xu = 0,
∂tu+ u∂xu+ 1
ρ∂xp = 0,
∂t(ρe) + ∂x(ρeu) + p∂xu = 0,
and obtain easily two basic transport equations for the motion ofρ and ρe as
∂tρ + u∂xρ = 0,
∂t(ρe) + u∂x(ρe) = 0.
By inserting the equation of state (6) into the latter one, we find an alternative description of the energy equation
∂t
µp+ pH
0 + ρeH
¶ + u∂x
µp+ pH
0 + ρeH
¶
= 0. (9)
To see how the pressure would retain in equilibrium as it should be for this model problem, we expand (9) into the form
1
0(∂tp+u∂xp)+ p
·
∂t
µ1 0
¶ +u∂x
µ1 0
¶¸
+
·
∂t
µpH
0 +ρeH
¶ +u∂x
µpH
0 +ρeH
¶¸
= 0.
Applying the assumed state of the pressure equilibrium to the above equation comes to a simplier expression as
p
·
∂t
µ1 0
¶ + u∂x
µ1 0
¶¸
+
·
∂t
µpH
0 + ρeH
¶ + u∂x
µpH
0 + ρeH
¶¸
= 0.
Now since this equation should hold for any p in the physical space, it implies right away that the terms in the square bracket of the equation should be vanished simultaneously, yielding a system of the following two equations:
∂t
µ1 0
¶ + u∂x
µ1 0
¶
= 0, (10a)
∂t
µpH
0 + ρeH
¶ + u∂x
µpH
0 + ρeH
¶
= 0. (10b)
It is important to note that in order to have the correct pressure equilibrium in (9) near the interfaces, (10a) and (10b) are the two key equations that should be satisfied for any given expression of0, pH, and eHin the Mie–Gr¨uneisen equation of state (6). From them, for a class of real materials modeled analytically by (6) at least, we may continue to work out suitable conditions for the further details of the related material parameters.
Consider the current two-phase flow application with0, pH, and eHdefined by (7), for example. Equations (10a) and (10b) now have the form
∂t
µ1− bρ γ − 1
¶ + u∂x
µ1− bρ γ − 1
¶
= 0, (11a)
∂t
µγ − bρ
γ − 1B + 2− γ − bρ γ − 1 aρ2
¶ + u∂x
µγ − bρ
γ − 1 B +2− γ − bρ γ − 1 aρ2
¶
= 0, (11b)
in a respective manner. It is obvious that, in addition to (11a) and (11b), we need to impose two supplementary conditions so as to have enough equations for the four material quantities:
γ , a, b, and B. In our approach, this is done quite easily by simply splitting (11a) into the following two parts,
∂t
µ 1 γ − 1
¶ + u∂x
µ 1 γ − 1
¶
= 0, (11c)
∂t
µ bρ γ − 1
¶ + u∂x
µ bρ γ − 1
¶
= 0, (11d)
and also (11b) into the terms
∂t
µγ − bρ γ − 1 B
¶ + u∂x
µγ − bρ γ − 1 B
¶
= 0, (11e)
∂t
µ2− γ − bρ γ − 1 aρ2
¶ + u∂x
µ2− γ − bρ γ − 1 aρ2
¶
= 0. (11f )
Having done so, we arrive at a primitive form of the transport equations (11c)–(11f) for the variables 1/(γ − 1), bρ/(γ − 1), B(γ − bρ)/(γ − 1), and aρ2(2 − γ − bρ)/(γ − 1). We note that with them for this interface only problem it is sufficient to have all the quantities γ , a, b, and B determined at all times, provided that the initial condition has been properly set for the computation; see Section 3.1.1.
Up to this point, our discussion stresses only on an approach that is capable of maintaining the pressure in equilibrium for a model interface only problem. Since in practice we are
interested in shock wave problems as well, we should thus take the equations, i.e., (11c)–
(11f ), in a form such thatγ , a, b, and B remain unchanged across both shocks and rarefaction waves. Concerning this, it is easy to see that withγ governed by (11c) there is no problem to do so (cf. [1, 63]). For b andB, however, due to the appearance of the linear factor of ρ in (11d) and (11e), it turns out that, in a time when such a scenario occurs, for consistency these two equations should be modified in such a way that each of them reduces to the basic mass conservation law of the fluid mixture. The derivation of the modification is simple.
Without going into the details, we write down the corrected version of the corresponding equation as follows:
∂t
µ bρ γ − 1
¶ + ∂x
µ bρ γ − 1u
¶
= 0, (12a)
∂t
µγ − bρ γ − 1 B
¶ + ∂x
µγ − bρ γ − 1 Bu
¶
= µ γ
γ − 1B
¶
∂xu. (12b)
We next come to the discussion for the quantity a described by (11f). In this situation, if we assume the proper smoothness of the nonlinear factors ofρ in the equation (such as in the case of rarefaction waves), as in (12a) and (12b) for the quantities b andB, respectively, we may derive an equation of the form
∂t
µ2− γ − bρ γ − 1 aρ2
¶ + ∂x
µ2− γ − bρ γ − 1 aρ2u
¶
= −
µ2− γ − 2bρ γ − 1 aρ2
¶
∂xu (12c)
that admits the desired constancy of a and the mass conservation as well. It is important to note that this is not the case, however, if the smoothness assumption onρ is violated; this happens in the event of shock waves, for there is no way to differentiate the discontinuous termsρ2andρ3yielding the partial derivatives onρ itself. Despite its apparent difficulty to fulfill the requirement of a across both shocks and rarefaction waves, (12c) is the correct form to be used in the model so that the fluid mixture aρ2(2 − γ − bρ)/(γ − 1) can be solved and applied readily to the computation of the pressure; see (14) below.
Here, rather than using (12c), we introduce a simplier linear advection equation of the form
∂ta+ u∂xa= 0 (12d)
for the motion of the mixture of a. Clearly, (12d) gives an accurate description of a in all the solution regions discussed above. In practice, it is a good one to use in the algorithm for numerical approximation, see Section 3.4 for some sample results.
To sum up, as in [63], we use the term effective equations to describe the set of equations that govern the motion of the material-dependent mixtures of the problems. In the current case, there consists of a couple sytem of five equations, i.e., (11c) and (12a)–(12d). It should be noted that when further taking the numerical aspect of the model approximation into consideration, with the full Euler equations, these effective equations are the proper ones to use for practical problems, but not the other form of the equations; see Section 3.3 and [60, 63] for the details.
Putting all the things together, with the modified van der Waals equation of state (5), the model equations we propose to solve one-dimensional multicomponent pro- blems are
∂tρ + ∂x(ρu) = 0
∂t(ρu) + ∂x(ρu2+ p) = 0
∂t(ρE) + ∂x(ρEu + pu) = 0
∂t
¡ bρ
γ − 1
¢+ ∂x
¡ bρ
γ − 1u¢
= 0
∂t
¡γ − bρ
γ − 1B¢ + ∂x
¡γ − bρ
γ − 1Bu¢
=¡ γ B
γ − 1
¢∂xu
∂t
¡2− γ − bρ
γ − 1 aρ2¢ + ∂x
¡2− γ − bρ
γ − 1 aρ2u¢
= −¡2− γ − 2bρ
γ − 1 aρ2¢
∂xu
∂t
¡ 1
γ − 1
¢+ u∂x
¡ 1
γ − 1
¢= 0
∂ta+ u∂xa= 0;
(13)
this gives us eight equations in total to be solved in one space dimension that is nicely independent of the number of fluid components involved in the problem. It is clear that in this system the first three are the Euler equations which are used to make certain the conservation of the basic fluid mixtures,ρ, ρu, and ρE, while the remaining are the effective equations that are introduced to ensure the correct mixing of the problem-dependent material variables near the interfaces. With a system expressed in this way, there is no problem to compute all the state variables of interest, including the pressure from the equation of state
p=
·
ρE −(ρu)2 2ρ −
µγ − bρ γ − 1B
¶
−
µ2− γ − bρ γ − 1 aρ2
¶¸Áµ 1
γ − 1− bρ γ − 1
¶ . (14)
Note that (13) reduces to a well-testedγ -based model for a stiffened gas when a = b = 0 (cf. [60, 63, 64]) and for a polytropic gas whenB = 0 as well (cf. [1]). As before, the proposed system (13) is not written in the full conservation form, but is rather a quasi-conservative system of equations. In addition, the nonzero terms on the right-hand side of the fifth and sixth equations of (13) should be viewed as an integrated part of the whole system, but not be considered as a source term. This fact can be realized easily by formulating the model as a quasi-linear system of equations (of course, with the assumption of a proper smoothness of the solutions), namely,
∂tq+ A(q)∂xq= 0 (15a)
with the state vector
q=
·
ρ, ρu, ρE, bρ
γ − 1,γ − bρ
γ − 1 B,2− γ − bρ
γ − 1 aρ2, 1 γ − 1, a
¸T
(15b)
and the matrix
A(q) =
0 1 0 0 0 0 0 0
K− u2 u(2 − 0) 0 p0 −0 −0 −p0 0
u(K − H) H − u20 u(0 + 1) up0 −u0 −u0 −up0 0
−ϕu ϕ 0 u 0 0 0 0
ϕuB −ϕB 0 0 u 0 0 0
−χu χ 0 0 0 u 0 0
0 0 0 0 0 0 u 0
0 0 0 0 0 0 0 u
. (15c)
Here K= 0u2/2, H = E + (p/ρ), ϕ = b/(γ − 1), and χ = aρ(4 − 2γ − 3bρ)/(γ − 1).
Recall that0 is the Gr¨uneisen coefficient, see (7), which gives different values for different fluids.
It is easy to show that for each physically relevant value of the state variables q defined in the region of thermodynamic stability (see Section 2), the eigenstructure of the matrix A is possessed of real eigenvalues
3 = diag(λ1, λ2, . . . , λ8) = diag(u − c, u, u + c, u, . . . , u) (16a) and a complete set of eigenvectors of the form
R= (r1, r2, . . . , r8) =
1 1 1 0 0 0 0 0
u− c u u+ c 0 0 0 0 0
H− uc 12u2 H+ uc −p 1 1 p 0
ϕ 0 ϕ 1 0 0 0 0
ϕB 0 ϕB 0 1 0 0 0
χ 0 χ 0 0 1 0 0
0 0 0 0 0 0 1 0
0 0 0 0 0 0 0 1
(16b)
with Ark= λkrk. Thus, (15a) is a hyperbolic system of partial differential equations and so is our multicomponent model (13). Regarding discontinuous solutions of the system, such as shock waves or contact discontinuities, (13) has the usual form of the Rankine–Hugoniot jump conditions across the waves; see Section 3.2.1 for more details.
With these comments, it should be sensible to use the proposed model for practical computations. The numerical method to be discussed in Section 3.3 is a consistent approxi- mation of the model that gives excellent results for a wide variety of problems as illustrated in Section 3.4.
3.1.1. Initialization of fluid-mixture cells. Consider a typical multicomponent setting in which there are m different fluids in a grid cell, and each of them occupies a distinct region with a volume-fraction function Y(i)in relation to it, for i= 1, 2, . . . , m. Here by the standard assumption, we have Y(i)∈ [0, 1] andPm
i=1Y(i)= 1. Suppose that for each component i the state variables such asρ(i), u(i), p(i), γ(i), a(i), b(i), andB(i)are known a priori. The objective is to give a proper definition of the fluid mixtures so that they can be used as an initial condition for our model equations (13) to the computations.
In the algorithm, we follow a common practice by evaluating the mixing statesρ, ρu, and ρe as a volume-weighted sum over the set of components ρ(i), ρ(i)u(i), andρ(i)e(i)given above,
ρ ρu ρe
=Xm
i=1
Y(i)
ρ(i) ρ(i)u(i) ρ(i)e(i)
. (17)
With this result, we then compute the mixture of the total energy byρE = ρe + (ρu)2/(2ρ);
this completes the definition of the conservative variables for the first three equations of the model.
To find the initial fluid mixtures, 1/(γ − 1), bρ/(γ − 1), B(γ − bρ)/(γ − 1), and aρ2(2 − γ − bρ)/(γ − 1), for the next four equations, we use the equation of state (5), where written as a function of the volume fraction it reads
µ1− bρ γ − 1
¶
(p + B + aρ2) + B − aρ2= ρe = Xm
i=1
Y(i)ρ(i)e(i)
= Xm
i=1
Y(i)
·µ1− b(i)ρ(i) γ(i)− 1
¶
(p(i)+ B(i)+ a(i)(ρ(i))2) + B(i)− a(i)(ρ(i))2
¸ .
By taking a similar approach as employed in Section 3.1 for the derivation of the effective equations, it comes out easily a splitting of the above equation into the form
γ − 11 bρ γ − 1 γ − bρ
γ − 1B
2− γ − bρ γ − 1 aρ2
= Xm
i=1
Y(i)
γ(i)1− 1
b(i)ρ(i) γ(i)− 1 γ(i)− b(i)ρ(i)
γ(i)− 1 B(i)
2− γ(i)− b(i)ρ(i)
γ(i)− 1 a(i)(ρ(i))2
, (18)
where in the process of splitting the terms the pressure p is chosen to satisfy the relation as µ1− bρ
γ − 1
¶ p=
Xm i=1
Y(i)
·µ1− b(i)ρ(i) γ(i)− 1
¶ p(i)
¸
(19)
(cf. [63] also for a simplier case of how this is done). With (18), it is easy to see that when each of the partial pressures is in equilibrium within a grid cell, the pressure p acquired from (19) would remain in equilibrium also; i.e., p= p(i), for i= 1, 2, . . . , m. Furthermore, from (18), we are able to obtain an explicit expression of the material-dependent parameters γ, a, b, and B in terms of the volume-fraction function and the original set of data. The results are
γ = 1 + 1 ,Ã m
X
i=1
Y(i) γ(i)− 1
!
, (20a)
b=
" m X
i=1
Y(i) Ã
b(i)ρ(i) γ(i)− 1
!#,"Ã m X
i=1
Y(i)ρ(i)
!Ã m X
i=1
Y(i) γ(i)− 1
!#
, (20b)
B =
" m X
i=1
Y(i)
Ãγ(i)− b(i)ρ(i) γ(i)− 1 B(i)
!#,"
1+ Xm
i=1
Y(i) Ã
1− b(i)ρ(i) γ(i)− 1
!#
, (20c)
a=
" m X
i=1
Y(i) Ã
2− γ(i)− b(i)ρ(i)
γ(i)− 1 a(i)(b(i))2
!# ,( Ã m X
i=1
Y(i)ρ(i)
!2
×
"
−1 + Xm
i=1
Y(i) Ã
1− b(i)ρ(i) γ(i)− 1
!#)
. (20d)
We note that as in the continuous counterpart (12c), (20d) is not used in practice for setting the initial mixture of a, but is done by the formula a=Pm
i=1Y(i)a(i)instead.
To end this subsection, it should be mentioned that by following the same approach introduced in [63] we may reformulate the model (13) into its variant form, the so-called volume-fraction model, that is robust when both the set of governing equations and the type of equations of state are different from one fluid component to the others, separated by the interfaces. To keep the presentation simple and clear, we omit the discussion of that model here, but refer the reader to [64, 65] for the details.
3.2. Approximate Riemann Solvers
Before describing numerical methods to solve (13), we pause to discuss the construction of the Riemann problem solutions which is one of the major steps in our multicomponent algorithm. For comparison purposes, we present two popular approaches for the resolution of the Riemann problem with piecewise constant data qLand qR to the left and right of the interface.
3.2.1. Shock-only solver. First, we are concerned with a shock-only 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. [3, 10, 13]). In this approach, the key step is to find the midstate (u∗, p∗) 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 p∗:
h(p∗) = u∗R(p∗) − u∗L(p∗) = 0. (21) Here u∗Land u∗Rare the velocities defined by connecting the states along the 1-shock and 3-shock curves, respectively,
u∗L(p) = uL− p− pL
ML(p), u∗R(p) = uR+ p− pR
MR(p),
with Mι denoting the Lagrangian shock speed, forι = L or R. In the current application with the modified van der Waals equation of state (5), when a= 0 (i.e., the vanishing of the aρ2term), we may compute Mιdirectly by evaluating the formula
Mι2(p) = Cι2
· 1+
µγι+ 1 2γι
¶µp+ Bι
pι+ Bι − 1
¶ ¸ ,
where Cι= ριcιis the Lagrangian sound speed (cf. [73]). Note that this is as a result derived from the Rankine–Hugoniot jump condition across the shock waves,
Mι2(e∗ι− eι) =¡
p2∗− p2ι¢±
2, (22)
with e∗ι= e(p∗, ρ∗ι) a function of the midstate density ρ∗ι, a quantity that is related to p, pι, and Mιin the following way
ρ∗ι(p) =
·
ρι−1− p− pι
Mι2(p)
¸−1 ,
while in the more general case when a6= 0, there is not such a close form solution available for Wι. Instead, we need to solve (22) iteratively for Mι; this is a typical thing to do when using the method for real gases (cf. [13]).
When applying a standard root-finding approach such as the secant method to (21), we have a 2-step iteration scheme as follows,
p(n+1)∗ = p∗(n)− ¯¯p∗(n)− p(n−1)∗ ¯¯
¯¯u(n)∗L − u(n−1)∗L ¯¯ + ¯¯u(n)∗R− u(n−1)∗R ¯¯£
u(n)∗R− u(n)∗L¤
, (23)
where u(n)∗ι = u∗ι[ p(n)∗ ], for ι = L or R, and n = 1, 2, . . . (until convergence). With a suitable choice of the starting values p∗(0) and p∗(1), method (23) typically converges to the exact solution p∗at a superlinear rate [31]. For gas dynamics, it is a common practice to set p∗(0) and p∗(1)by
p∗(0)= pRCL+ pLCR− (uR− uL)CLCR
CL+ CR
, u(0)∗L = uL− p(0)∗ − pL
CL
, u(0)∗R= uR+ p(0)∗ − pR
CR
, (24)
p∗(1)= pRM(0)L + pLM(0)R − (uR− uL)M(0)L MR(0) M(0)L + M(0)R
,
with Mι(0)= Mι[ p∗(0)]. After a satisfactory convergence of the scheme, we may then calculate u∗by
u∗= pL− pR+ uLML(p∗) + uRMR(p∗) ML(p∗) + MR(p∗) .
We note that alternatively we may use a 1-step Newton method for the solution of (21).
Since the derivation of the scheme is more involved due to the need to compute the derivative term d p∗/du∗, we do not discuss the method here, but refer the reader to [13] for more details.
Figure 2 shows a typical solution structure of the 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
λ1= u∗− ML(p∗)
ρ∗L(p∗), λ2= u∗, λ3= u∗+ MR(p∗)
ρ∗R(p∗), (25a)
FIG. 2. Typical solution structure of the Riemann problem for our multicomponent model (13). Note that in a shock-only approximation of the approximate Riemann solver, the rarefaction wave is replaced by an entropy- violating shock.
with the jumps across each of them computed by the difference between the states to the left and right of the discontinuity,
W1 = q∗L− qL, W2= q∗R− q∗L, W3= qR− q∗R. (25b) Wave propagation methods are based on using these propagating discontinuities to update the cell averages in the cells neighboring each interface.
It is true that no matter what iterative method is employed to the solution of (21), the approach is quite expansive as compared to the approximate solver of Roe described next.
Nevertheless, there are various situations where this approach is worthwhile and can provide more accurate results than the Roe solver does. This includes some examples shown in [18, 41, 56] and many difficult problems with strong shock waves and stiff equations of state. Moreover, it is a straightforward matter to generalize the approach that covers the case with surface tension effect across the interfaces; results on this aspect will be reported elsewhere.
3.2.2. Roe solver. In a Roe’s approximate Riemann solver, we replace the nonlinear system (13) with data qLand qRby a linear system of the form
∂tq+ ˆA(qL, qR)∂xq = 0. (26) Here ˆA(qL, qR) is a constant matrix that depends on the initial data and is a local linearization of the matrix A in (15c) about an average state. To find that matrix, as it is often done in many other Roe solvers (cf. [9, 22, 23]), we want to seek an average state such that the difference of the fluxes in the conservation part of (13) (i.e., the first four equations of the system) are equal to the respective first order approximations of the flux differences. That is, 1F(i)= (FR− FL)(i)= [ ˆA(qL, qR)(qR− qL)](i)= [ ˆA(qL, qR)1q](i), (27) for i = 1, 2, 3, 4, where F ∈ R4is the usual definition of the fluxes for conservation laws, and1F(i)is the i th component of1F. With that, it is a straightforward matter to obtain the results for ˆu, ˆH, and ˆϕ by the standard “Roe-averaging” approach; i.e., for a given pair (ρL, ρR), the average state for a quantity z is defined by
ˆz=√ρLzL+ √ρRzR
√ρL+ √ρR
. (28)
Note that in the process of the derivation, as in [63] we have chosen the averages(d1/0) and(dp/0) based on (28) so that the expression
1p =· dµ 1 0
¶ 1
µp 0
¶
−µd p 0
¶ 1
µ1 0
¶¸Á dµ 1 0
¶2
is satisfied approximately. With that we set ˆp= d(p/0)/ d(1/0) and ˆ0 = 1/ d(1/0). To finish the construction of ˆA(qL, qR), we still need to find the averages of B and χ. Since there is no unique way to do so, we might as well compute ˆB and ˆa using (28) and set ˆχ = [(2/ ˆ0) − 2 − ˆϕ ˆρ]ˆa ˆρ, where ˆρ = √ρLρR. Numerical results shown in Section 3.4 indicate that the set of average states described here is a good one to use for practical multicomponent problems.
The solution of the linear problem (26) consists of eight discontinuities propagating at constant speeds (for a system of eight equations). The jump across each discontinuity is a multiple of the eigenvector of the matrix ˆA, and the propagating speed is the corresponding eigenvalue. We thus have
1q = qR− qL= X8 k=1
αˆkˆrk, (29)
where ˆrkis the kth eigenvector of ˆA; see (16a) and (16b). The scalar ˆαkgives the strength across the discontinuity that can be determined easily from (29). We find
αˆ2= 1q(1)+ 0ˆ ˆ c2
·
−ˆu2
2 1q(1)+ ˆu1q(2)− 1q(3)+ ˆp¡
1q(7)− 1q(4)¢
+ 1q(5)+ 1q(6)
¸ , αˆ3= 1
2 ˆc
£(ˆc − ˆu)1q(1)+ 1q(2)− ˆc ˆα2
¤,
αˆ1= 1q(1)− ˆα2− ˆα3, αˆ4= 1q(4)− ˆϕ(1ρ − ˆα2), (30) αˆ5= 1q(5)− ˆϕ ˆB(1ρ − ˆα2), ˆα6= 1q(6)− ˆχ(1ρ − ˆα2),
αˆ7= 1q(7), αˆ8= 1q(8), where ˆc=p
0[ ˆH − (ˆuˆ 2/2) + ˆp ˆϕ − ˆϕ ˆB − ˆχ] is the speed of sound.
Note that in this Riemann solution, except the discontinuities for ˆλ1= ˆu − ˆc, and ˆλ3= ˆu + ˆc, all the other discontinuities (six of them) are propagating at the same speed ˆu. For practical purposes, we may view these discontinuities as a single one with the operator W2 defined by combining all the jumps across the ˆλ2 wave family; i.e., set W2= ˆα2ˆr2+P8
k=4αˆkˆrk. Clearly, doing so removes the effect of the wave families ˆλ4
to ˆλ8from the solution. With this notation, we also writeWk= ˆαkˆrkto represent the jump across the k-wave for k= 1 or 3.
3.3. Wave Propagation Methods
We use the high-resolution method based on a wave-propagation viewpoint to compute approximate solutions of our multicomponent model introduced in Section 3.1. The method is a variant of the fluctuation-and-signal scheme of Roe [58, 59] in that we solve the 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 (cf. [34, 39]).
For simlicity, we describe the method on a uniform grid with fixed mesh spacing1x, but the method can be extended quite easily to a nonuniform and time-varying grid as well (cf.
[35, 40]). We use a standard finite-volume formulation in which the value Qnjapproximates the cell average of the solution over the grid cell [xj, xj+1] at time tn:
Qnj ≈ 1 1x
Z xj+1 xj
q(x, tn) dx.
The time step from the current time tnto the next tn+1is denoted by1t.
3.3.1. First order method. In this setup, a first order accurate version of the method in wave-propagation form is a Godunov-type scheme that can be written as
Qnj+1= Qnj− 1t 1x
mw
X
k=1
(λ−kWk)nj+1+ (λ+kWk)nj, (31)
whereλk∈ R and Wk∈ Rmare solutions of the kth wave family, for k= 1, 2, . . . , mw, ob- tained from solving the Riemann problems at cell interfaces xjand xj+1, λ− = min(λ, 0), and λ+ = max(λ, 0). Clearly, the method belongs to a class of upwind schemes (cf.
[23, 36]), and it will be shown next that the method is quasi-conservative in the sense that when applying the method to (13) not only the conservation equations but also the transport equations are approximated in a consistent manner by the method with the chosen Riemann solver.
To demonstrate that, we first analyze our method for an interface only problem as de- scribed in Section 3.1. Without loss of generality, we consider a single Riemann problem where at cell interface xjthe initial data consists of uniform pressure p0and constant par- ticle velocity u0 to the left and right of the interface, but with jumps on the other state variables of q. Assuming a positive velocity u0 > 0, for example. If the problem is solved by using a shock-only Riemann solver (see Section 3.2.1), from (31) the cell average Qnj would be updated by
Qnj+1= Qnj− 1t
1x(λ2W2)nj, (32a)
or equivalently by
ρ ρu ρE
bρ γ − 1 γ − bρ
γ − 1B
2− γ − bρ γ − 1 aρ2
1 γ − 1
a
n+1
j
=
ρ ρu0
ρE
bρ γ − 1 γ − bρ
γ − 1B
2− γ − bρ γ − 1 aρ2
1 γ − 1
a
n
j
− 1t 1xu0
1ρ u01ρ 1ρE 1γ − 1bρ 1γ − bργ − 1B 12− γ − bργ − 1 aρ2
1γ − 11 1a
n
j
, (32b)
when expressing (32a) in terms of the solution states of the problem. Noting that in this case the difference operator 1 is simply applied to the Riemann data Qnj−1 and Qnj on
