An Efficient Shock-Capturing Algorithm for Compressible Multicomponent Problems

(1)

JOURNAL OF COMPUTATIONAL PHYSICS142, 208–242 (1998) ARTICLE NO. CP985930

An Efficient Shock-Capturing Algorithm for Compressible Multicomponent Problems

Keh-Ming Shyue

Department of Mathematics, National Taiwan University, Taipei, Taiwan, Republic of China E-mail: shyue@math.ntu.edu.tw

Received January 6, 1997; revised January 22, 1998

A simple shock-capturing approach to multicomponent flow problems is developed for the compressible Euler equations with a stiffened gas equation of state in multiple space dimensions. The algorithm uses a quasi-conservative formulation of the equations that is derived to ensure the correct fluid mixing when approximating the equations numerically with interfaces. Aγ -based model and a volume-fraction model have been described, and both of them are solved using the standard high-resolution wave propagation method for general hyperbolic systems of partial differential equations. Several calculations are presented with a Roe approximate Riemann solver that show accurate results obtained using the method without any spurious oscillations in the pressure near the interfaces. Convergence of the computed solutions to the correct weak ones has been verified for a two-dimensional Richtmyer–Meshkov unstable interface problem where we have performed a mesh-refinement study and also shown front-tracking results for comparison. ° 1998 Academic Pressc

1. INTRODUCTION

Our goal is to present a simple approach to multicomponent flow of general compressible materials in more than one dimension. We use the Euler equations of gas dynamics as a model system, and consider problems with the so-called “stiffened” gas equation of state for approximating materials including compressible liquids and solids [29, 50]. The algorithm uses a quasi-conservative formulation of the equations proposed by Abgrall [1] to ensure a consistent approximation of the energy equation near the interfaces where regions of different fluid components are separated. Aγ -based model is therefore derived that extends the work of Abgrall [1] from polytropic gases in one dimension to a stiffened gas, and also to multiple dimensions. We give a new formulation of the resulting model in expression of the volume fraction that is more robust for two-component flow problems. This will be discussed further in Sections 2 and 5.

208 0021-9991/98 $25.00

Copyright c° 1998 by Academic Press All rights of reproduction in any form reserved.

(2)

MULTICOMPONENT ALGORITHM FOR COMPRESSIBLE FLOW 209

We use the high-resolution wave propagation method developed by LeVeque [40, 41, 43]

to solve the proposed multicomponent models. This method is a variant of the fluctuation- and-signal method of Roe [59, 60], and has been widely used in many applications including the single-component fluid of ideal gases. The main idea behind the method has recently been implemented in the software package CLAWPACK (Conservation LAWs PACKage) as the underlying integration routine [46, 47]. The current use of the method is just an easy extension of the previous one from single-component to multicomponent problems. It is an efficient and yet accurate scheme without any spurious oscillations in the pressure near an interface as illustrated by numerical results presented in this paper.

We will only briefly review and describe the method in a shock-capturing framework in one dimension, see Section 3. Extensions of the method to front-tracking and to two dimensions are straightforward also by following the procedures outlined in [44, 45], for example. We will not discuss a front-tracking version of the method here, but fo- cus on the more fundamental shock-capturing algorithm and validate its use via numeri- cal experimentations. Some preliminary results obtained using the front tracking method may be found in [65, 66] for two-dimensional unstable fluid interface problems such as Rayleigh–Taylor and Richtmyer–Meshkov instabilities, see Section 6 for an example also.

Generalization of the approach to three dimensions can be made in a similar manner, but requires more programming work, especially in regard to a front-tracking method [25].

Clearly for real applications the use of a stiffened gas equation of state that appears in an analytical formula (i.e., the constitutive relation (2) in Section 2) represents only a limited number of materials of practical importance [29, 50]. However, there are some problems of sufficient interest and difficulty that the development of a multicomponent algorithm for this equation of state is worthwhile, particularly since in some cases it is relatively easy to compute the exact solutions and check accuracy of the method.

Numerous numerical methods have been developed over the years to handle multicomponent flow problems. Consider a non-reacting ideal gas flow, for example. One popular approach among them is to solve an extended system of equations in which additional conservation equations are introduced to the original Euler equations to describe the conservation of each fluid component separately. Methods of this type, in particular a shock- capturing version of the method, often fail to maintain pressure equilibrium for grid cells near interfaces where two or more fluid components are mixed. Some representative examples that exhibit this erroneous phenomenon are given in [34] for the use of aγ -based (see Fig. 1 also) and a level-set model, and in [9, 12, 71] for the use of a mass-fraction model.

The exception is the method explained by Jenny et al. [33] that the fluxes obtained using conservative Euler solvers are modified in a suitable way to avoid the occurrence of the pressure errors generated near the interfaces. It is unclear how to extend the method to a situation other than the ideal gas flows, however.

Another approach introduced by Karni [35, 36] is to solve the Euler equations separately on each side of the interface using a method designed for a single-component flow, while the interface is dealt with in a different manner using a pressure evolution equation derived from the energy equation. Despite the fact that the method is not exactly conservative at the interface, reasonable results are obtained using this approach in conjunction with either standard level-set or mass-fraction formulation of ideal gases. Extension of the method to a thermally perfect gas was done recently by Fedkiw et al. [21] in one dimension. It should be noted that Cocchi et al. [12] devised a rather similar method of this kind that employs a

(3)

210 KEH-MING SHYUE

linear interpolation technique near the interface instead; some one-dimensional results are shown for a stiffened gas flow.

Our method to modeling multicomponent flow of general compressible materials is mo- tivated by the work of Abgrall [1] in that, based on the physical principles and from the energy equation, we derive the effective equations (see Section 2) for the mixture of material- dependent quantities near the interface. In this method, with the stiffened gas equation of state, we take these equations to be of the form that do not vary their solutions across the shock and rarefaction waves as well. Combining the resulting effective equations (i.e., Eqs. (9)) to the Euler equations yields a model system that is not written in the full conservation form, but is rather a quasi-conservative system of equations. Abgrall [1] solves a system of this kind using a predictor-corrector method, while we use the high-resolution wave propagation method that gives an efficient implementation of the algorithm. In prin- ciple, when properly modified, it is possible to employ the state-of-the-art shock-capturing methods for hyperbolic systems of conservation laws for solving the model equations as well; see [67] for an example that generalizes the MUSCL scheme [13, 72] to this specific application and also to a van der Waals gas.

There are many other multicomponent approaches available in the literature. Some typical ones are the level-set methods [19, 52], front-tracking methods [25, 28], volume-of-fluid methods [14, 51, 55, 71], and the BGK-based method [74].

This paper is organized as follows. In Section 2, we begin by discussing the basic computational models in one dimension that govern the motion of materials characterized by a stiffened gas equation of state in a multicomponent environment. We will give two different formulations of the model equations that are both applicable for practical computations of multicomponent problems. In Section 3, we briefly review an approximate Riemann solver of Roe, and analyze a first order numerical methods based on a wave-propagation approach with an application to our multicomponent models proposed in Section 2. The results of some one-dimensional tests are given in Section 4 that validate this approach. The basic idea of the algorithm is then extended to multiple dimensions in Section 5, and some two-dimensional numerical results are presented. Section 6 shows results for a Richtmyer–

Meshkov unstable interface problem where a mesh-refinement study is performed to check convergence of the computed solutions.

2. DERIVATION OF MODEL EQUATIONS

In one dimension, the single-component Euler equations of gas dynamics take the form

∂

∂t



 ρ ρu ρE



 + ∂∂x



 ρu

ρu²+ p (ρE + p)u



 = 0, (1)

whereρ is the density, u is the velocity, p is the pressure, and E is the total energy per unit mass. We assume a compressible material that the internal energy per unit mass, denoted by e, satisfies the “stiffened” gas equation of state,

ρe = p+ γ p_∞

γ − 1 , (2)

(4)

and E= e + u²/2. Here γ is the usual ratio of specific heats (γ > 1), and p∞ is a pre- scribed pressure-like constant; these values can be used to describe the material property of interests and can be determined from laboratory experiments via an empirical fit [31, 49].

For example, for water we haveγ ∼= 5.5, p∞= 4921.15 bars [12], and for tungsten we getγ ∼= 3.14, p∞= 1.0 Mbar [54]. Note a stiffened gas reduces to a polytropic gas when p_∞= 0. The three components of Eqs. (1) express the conservation of mass, momentum, and energy, respectively [18].

We are interested in the simulation of multicomponent flow problems. For the equations, we take a popular approach by considering the Euler Eqs. (1) as a model system, see [33, 74] for the use of other governing equations. Our goal here is to derive computational models that may prevent pressure oscillations near the interfaces, when solving the problem numerically with standard shock-capturing methods.

2.1. Preliminary. To begin, suppose that there are m different fluid components in a grid cell, and each of them occupies a distinct region with a volume-fraction function Y⁽ⁱ⁾in relation to it, for i= 1, 2, . . . , m. Here by the standard assumptions we have Y⁽ⁱ⁾

∈ [0, 1] andPm

i=1Y⁽ⁱ⁾= 1. Suppose that for each component i the state variables such as ρ⁽ⁱ⁾, u⁽ⁱ⁾, p⁽ⁱ⁾, γ⁽ⁱ⁾, and p⁽ⁱ⁾_∞are known a priori. The objective is to define the mixture of the pressure p as well as the conserved variablesρ, ρu, and ρE in a consistent manner within the cell for the Euler Eqs. (1). Note this step is necessary when we initialize the data for computations.

To accomplish this, we follow a common practice by settingρ, ρu, and ρE as a volume- weighted sum over the set of componentsρ⁽ⁱ⁾, ρ⁽ⁱ⁾u⁽ⁱ⁾, andρ⁽ⁱ⁾E⁽ⁱ⁾, for each separately,

ρ = Xm

i=1

Y⁽ⁱ⁾ρ⁽ⁱ⁾, ρu = Xm

i=1

Y⁽ⁱ⁾ρ⁽ⁱ⁾u⁽ⁱ⁾, ρE = Xm

i=1

Y⁽ⁱ⁾ρ⁽ⁱ⁾e⁽ⁱ⁾+1

2Y⁽ⁱ⁾ρ⁽ⁱ⁾u², (3) where u is the velocity mixture that can be computed easily by

u =ρu ρ =

Pm

i=1Y⁽ⁱ⁾ρ⁽ⁱ⁾u⁽ⁱ⁾ Pm

i=1Y⁽ⁱ⁾ρ⁽ⁱ⁾ .

To find the pressure mixture p, we need to use the equation of state (2). From (3), after some simple algebra, we find

p+ γ p∞

γ − 1 = ρe = Xm

i=1

Y⁽ⁱ⁾ρ⁽ⁱ⁾e⁽ⁱ⁾= Xm

i=1

Y⁽ⁱ⁾ Ã

p⁽ⁱ⁾+ γ⁽ⁱ⁾p_∞⁽ⁱ⁾ γ⁽ⁱ⁾− 1

!

. (4)

Note this gives us one equation that is not only for the mixture of p but also forγ and p_∞. Obviously, we need to choose two supplementary conditions so as to have p,γ , and p∞

defined well, leading to the full agreement of (4).

The basic idea of our approach is quite simple. We first split (4) into two parts by setting the terms

p γ − 1 =

Xm i=1

Y⁽ⁱ⁾p⁽ⁱ⁾

γ⁽ⁱ⁾− 1 and γ p_∞ γ − 1 =

Xm i=1

Y⁽ⁱ⁾γ⁽ⁱ⁾p⁽ⁱ⁾_∞

γ⁽ⁱ⁾− 1 ; (5)

(5)

212 KEH-MING SHYUE

this provides us with one condition right away. We then impose the following condition to the computation ofγ ,

1 γ − 1=

Xm i=1

Y⁽ⁱ⁾

γ⁽ⁱ⁾− 1. (6a)

Clearly whenγ is known, from (5), it is an easy matter to determine the unknowns p and p_∞. The results are

p= Ã _m

X

i=1

Y⁽ⁱ⁾p⁽ⁱ⁾ γ⁽ⁱ⁾− 1

!,Ã _m X

i=1

Y⁽ⁱ⁾ γ⁽ⁱ⁾− 1

!

(6b)

and

p_∞= Ã _m

X

i=1

Y⁽ⁱ⁾γ⁽ⁱ⁾p⁽ⁱ⁾_∞ γ⁽ⁱ⁾− 1

!,Ã 1+

Xm i=1

Y⁽ⁱ⁾ γ⁽ⁱ⁾− 1

!

. (6c)

Notice that in case each of the partial pressures p⁽ⁱ⁾ is in equilibrium within a grid cell the pressure p acquired from (6b) would remain in equilibrium also, i.e., p= p⁽ⁱ⁾, for i= 1, 2, . . . , m. This nice property on the mixture of p is in fact the main reason why we use (6a), but not other ways to computeγ , see [14, 55] for a more rigorous derivation. In addition, we get the total pressure p=Pm

i=1Y⁽ⁱ⁾p⁽ⁱ⁾when the grid cell is composed of the same set ofγ⁽ⁱ⁾but with different p⁽ⁱ⁾. We use (6c) for the mixture of p_∞in order to make sure that (4) is handled in a consistent manner.

It is worthwhile to mention that in case Y⁽ⁱ⁾represents a mass-fraction function of the i th component withρ⁽ⁱ⁾= ρY⁽ⁱ⁾, we would have Eq. (4) replaced by

p+ γ p_∞

γ − 1 = ρe = Xm

i=1

ρ⁽ⁱ⁾e⁽ⁱ⁾= Xm

i=1

Ã

p⁽ⁱ⁾+ γ⁽ⁱ⁾p⁽ⁱ⁾_∞ γ⁽ⁱ⁾− 1

! .

In this instance, it is a standard approach that setsγ according to (6a) (cf. [26, 39]).

Analogously, by following the same procedures as for the volume-fraction case, we find the results of p and p_∞,

p = Ã _m

X

i=1

p⁽ⁱ⁾ γ⁽ⁱ⁾− 1

!,Ã _m X

i=1

Y⁽ⁱ⁾ γ⁽ⁱ⁾− 1

! ,

p_∞= Ã _m

X

i=1

γ⁽ⁱ⁾p⁽ⁱ⁾_∞ γ⁽ⁱ⁾− 1

!,Ã 1+

Xm i=1

Y⁽ⁱ⁾ γ⁽ⁱ⁾− 1

! .

2.2. γ -based model. It is clear that the motion of the fluid mixtures, obtained from above, is governed by the Euler Eqs. (1). In the development of our multicomponent model, it is of great importance to first consider the case of an interface only problem where both of the pressure p and velocity u are constants in the domain, while the other variables such asρ, γ , and P_∞are having jumps across some interfaces.

(6)

To do so, we start out to write (1) in the following non-conservative form,

∂ρ

∂t + u∂ρ

∂x + ρ∂u

∂x = 0,

∂u

∂t + u∂u

∂x +1 ρ

∂p

∂x = 0,

∂

∂t(ρe) + ∂

∂x(ρeu) + p∂u

∂x = 0, and obtain easily equations describing the motion of the interfaces as

∂ρ

∂t + u∂ρ

∂x = 0, (7a)

∂

∂t(ρe) + u ∂

∂x(ρe) = 0. (7b)

With this, it is clear that for an interface the densityρ as well as the total internal energy ρe is evolved by the linear transport Eqs. (7a) and (7b).

To see how the pressure p would retain in equilibrium as it should be for this problem, we insert the equation of state (2) into (7b), and have

∂

∂t

µp+ γ p_∞ γ − 1

¶ + u ∂

∂x

µp+ γ p_∞ γ − 1

¶

= 0. (8)

By expanding (8), we may therefore write that equation as µ 1

γ − 1

¶·∂p

∂t + u∂p

∂x

¸ + p

·∂

∂t µ 1

γ − 1

¶ + u ∂

∂x µ 1

γ − 1

¶¸

+

·∂

∂t µγ p_∞

γ − 1

¶

+ u ∂

∂x µγ p_∞

γ − 1

¶¸

= 0.

Now the requirement that p be in equilibrium leads to the equation of the form

p

·∂

∂t µ 1

γ − 1

¶ + u ∂

∂x µ 1

γ − 1

¶¸

+

·∂

∂t µγ p_∞

γ − 1

¶ + u ∂

∂x µγ p_∞

γ − 1

¶¸

= 0.

Since this equation should hold for any p in the physical space, it implies that the terms in bracket of the above equation should be vanished simultaneously, yielding a system of two equations











∂

∂t µ 1

γ − 1

¶ + u ∂

∂x µ 1

γ − 1

¶

= 0

∂

∂t µγ p_∞

γ − 1

¶ + u ∂

∂x µγ p_∞

γ − 1

¶

= 0.

(9)

Note that these are the evolution equations that should be satisfied for the material-dependent variablesγ and p_∞in order to have the correct pressure behavior in (8) for the interface. For convenience, we call (9) the effective equations of the problem, where the initial condition of the equations are provided, for example, by (6a) and (6c) in a respective manner.

(7)

214 KEH-MING SHYUE

Of course, intuitively, there are many other ways that the effective equations can be rewritten, while still getting the correct pressure from (8) for this interface only problem.

One simple example among them is to write (9) as











∂

∂t µ 1

γ − 1

¶ + ∂

∂x µ 1

γ − 1u

¶

= 0

∂

∂t µγ p∞

γ − 1

¶ + ∂

∂x µγ p∞

γ − 1u

¶

= 0 ;

(10)

this is a legitimate one to use, for the velocity u is a constant for the problem (see Eqs. (25a)–(25c) also for other examples when p_∞= 0). But since in general we are interested in shock wave problems as well, with the stiffened gas equation of state (2), we should only take the effective equations in a form that do not vary their solutions across both of shocks and rarefaction waves. For this reason, it rules out immediately the use of (10) as the effective equations. When further taking the numerical aspect of the model equations into consideration, it turns out that with the full Euler Eqs. (1) the effective equations in a form of (9) are the proper ones to use for practical multicomponent problems; see Section 3 for the analysis of a numerical method that approximates the model equations, and Section 4 for numerical examples.

In summary, the model equations we propose to solve multicomponent problems with the stiffened gas equation of state are











∂ρ

∂t + ∂

∂x(ρu) = 0

∂

∂t(ρu) + ∂

∂x(ρu²+ p) = 0

∂

∂t(ρE) + ∂

∂x[(ρE + p)u] = 0

∂

∂t µ 1

γ − 1

¶ + u ∂

∂x µ 1

γ − 1

¶

= 0

∂

∂t µγ p∞

γ − 1

¶ + u ∂

∂x µγ p∞

γ − 1

¶

= 0.

(11)

Here in this system the first three equations are simply the Euler Eqs. (1) that are used to make certain the conservation ofρ, ρu, and ρE, while the last two equations are the effective Eqs. (9) of the problem, that are introduced to ensure the correct mixing behavior of the variablesγ and p∞on the interfaces.

Note the model Eqs. (11) is not written in the full conservation form, but is rather a quasi-conservative system of equations. For shock wave problems, this poses no problem at all, because for stiffened gasesγ and p_∞remain unchanged across genuinely nonlinear waves such as shocks or rarefactions, and hence we have the usual Riemann invariants and Rankine–Hugoniot jump conditions across the rarefaction and shock waves, respectively (cf. [26, 68]). At linearly degenerate waves such as interfaces where there may be jumps in γ and p∞, this again gives no problem since as we have seen in the previous discussion we would have the desired pressure equilibrium when the model is in use.

Of particular importance is the case of shock wave and interface interaction. It is known that, due to the nonlinearity of the problem, the pressure across the interface will be quite

(8)

different before and right after the wave interaction [7, 18]. Somewhat surprisingly, we find no major problem in using the model as well, see Section 4 for a representative numerical test. We emphasize that in this caseγ and p_∞are transported in a passive manner according to (9) along with the interface, and the variables such asρ, u, and p are dealt with in a conservative way from the conservation part of Eqs. (11) as usual.

With all these in mind, it should be sensible to use the proposed model as for practical computations. The numerical method to be described in Section 3 is a consistent approximation of the model that gives excellent results for a wide variety of problems as illustrated in Section 4. For convenience, we call this model aγ -based model to be distinct from the other one presented below. Note in [1] Abgrall used a slightly different starting point, but nevertheless he obtained the same set of model equations for polytropic gas problems when

p_∞= 0.

2.3. Volume-fraction model. It should be noted that we may reformulate the aboveγ - based model as a volume-fraction or a mass-fraction model that is also applicable to many multicomponent problems. To demonstrate the basic idea, we consider the case in using volume-fraction functions as an example.

Similar to the approach used in the construction of our γ -based model, we look for effective equations that may preserve the pressure equilibrium for an interface only problem.

Here with the volume-fraction notion of the states 1/(γ −1) and γ p∞/(γ −1) being defined by (6a) and (6c), the key step is to replace (9) by

∂

∂t Ã _m

X

i=1

Y⁽ⁱ⁾ γ⁽ⁱ⁾− 1

! + u ∂

∂x Ã _m

X

i=1

Y⁽ⁱ⁾ γ⁽ⁱ⁾− 1

!

= 0,

∂

∂t Ã _m

X

i=1

Y⁽ⁱ⁾γ⁽ⁱ⁾p_∞⁽ⁱ⁾ γ⁽ⁱ⁾− 1

! + u ∂

∂x Ã _m

X

i=1

Y⁽ⁱ⁾γ⁽ⁱ⁾p⁽ⁱ⁾_∞ γ⁽ⁱ⁾− 1

!

= 0.

After regrouping terms, we find the transport equation for each volume fraction Y⁽ⁱ⁾,

∂Y⁽ⁱ⁾

∂t + u∂Y⁽ⁱ⁾

∂x = 0 for i = 1, 2, . . . , m. (12) As before when the set of Y⁽ⁱ⁾is known, we may therefore computeγ and p_∞from (6a) and (6c). In effect, in a volume-fraction model, instead of using (9) we use (12) (m of them) as the effective equations of the problem. Note that in case Y⁽ⁱ⁾stands for the mass-fraction of the component i , we find the effective equations of a mass-fraction model that is of the same form as in (12). Because of the close connection between the two models, we devote our discussion to one of the models only, namely to the volume-fraction model.

For completeness, we write down the full set of equations for this volume-fraction model,











∂ρ

∂t + ∂

∂x(ρu) = 0

∂

∂t(ρu) + ∂

∂x(ρu²+ p) = 0

∂

∂t(ρE) + ∂

∂x[(ρE + p)u] = 0

∂Y⁽ⁱ⁾

∂t + u∂Y⁽ⁱ⁾

∂x = 0 for i = 1, 2, . . . , m;

(13)

(9)

216 KEH-MING SHYUE

this gives us totally m+ 3 equations to be solved. Since the derivation of the volume- fraction model comes closely out of theγ -based model, it can be shown that this model is as effective as theγ -based model. But for general multicomponent problems, the γ -based model is the preferred one to use, because the basic equations for the model stay as five, see (11), irrespective of the number of components involved in the problem.

3. NUMERICAL APPROXIMATION OF MODEL EQUATIONS

We use the high-resolution wave propagation method to compute solutions for our mul- ticomponent models introduced in Section 2. The method is a variant of the fluctuation- and-signal formulation of Roe [59, 60] in that we solve the Riemann problem at each cell interface, and use the resulting waves (i.e., discontinuities propagating at constant speeds) to update the solutions in neighboring grid cells (cf. [40, 41]).

For simplicity, 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 [44]. We use a standard finite-volume formulation that the value Uⁿ_j ∈ R^mapproximates the cell average of the solution over the grid cell [xj, xj+1] at time tn. The time step is denoted by1t.

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

Uⁿ_j⁺¹= Uⁿj − 1t 1x

Xm k=1

£(λ⁻kWk)ⁿj+1+ (λ⁺kWk)ⁿj

¤, (14)

whereλk∈ R and Wk∈ R^mare solutions obtained from solving Riemann problems at cell interfaces xj and xj+1, λ⁻ = min(λ, 0), and λ⁺ = max(λ, 0) (cf. [26]). It is easy to see that the method belongs to a class of upwind schemes, and as it will be shown below that the method is quasi-conservative in the sense that when applying the method to our multicomponent models not only the conservation equations but also the transport equations are approximated in a consistent manner by the method with the chosen Riemann solver.

Concerning stability, it is observed numerically that the method is stable and convergent under mesh refinement provided that the waves in the method affect only the cells adjacent to the interface during the time step; see Section 4 for numerical examples and also the results shown in [43].

3.1. Method withγ -based model. Consider the γ -based model as an example. We solve the Riemann problem at each cell interface x_jthat consists of (11) with piecewise constant data Uⁿ_j₋₁and Uⁿ_j on the left and on the right of the interface. Rather than computing the exact solution to this Riemann problem, which can be done by iterative procedures (cf. [12, 15, 54]) but is rather expensive, we use a generalized version of the approximate Riemann solver of Roe (cf. [58] and below) in most instances. This is much more efficient to compute than the exact Riemann solution, and provides a very accurate approximation of solution for smooth flows and also for moderate-strength shock waves. As long as the equation of state is not too stiff across the interfaces (which is the application considered here), the solution of the Roe Riemann solver gives a proper resolution of the contact discontinuity to be used in the method (14) for numerical approximation also. It is important to mention that in case we have a stringent set of data that involves strong shock waves and (or) stiff equation of states, we find it is advisable to use the exact Riemann solver so as to properly deal with

(10)

the nonlinear effect of the solution structure; see [67] for an example and also [20, 45, 56]

for other examples that the exact Riemann solver should be in use.

To implement Roe’s approximate Riemann solver, we first write (11) as a quasi-linear system of equations

∂q

∂t + A(q)∂q

∂x = 0 (15)

with

q=





 ρ ρu ρE

1 γ −1 γ p∞

γ −1







, A(q) =







0 1 0 0 0

³γ −3 2

´

u² (3 − γ )u γ − 1 χ 1− γ

³γ −1 2

´

u³− uH H − (γ − 1)u² γ u χu (1 − γ )u

0 0 0 u 0

0 0 0 0 u





 ,

where H = E + (p/ρ) is the enthalpy, and χ = −ι(γ − 1)²withι = p/(γ − 1). We then solve a linear problem

∂q

∂t + ˆA(qL, qR)∂q

∂x = 0 (16)

with initial data

q(x, 0) =

½qL for x left to the interface qR for x right to the interface,

where ˆA(qL, qR) is a constant matrix that depends on the initial data and is a local lineariza- tion of the matrix A about an average state. Here as it is often done in many other Roe solvers (cf. [10, 23, 26]), we want to seek an average state that the difference of the fluxes in the conservation part of Eqs. (11) are equal to the respective first order approximation of the flux differences. That is,

1F⁽ⁱ⁾= (FR− FL)⁽ⁱ⁾= [ ˆA(qL, qR)(qR− qL)]⁽ⁱ⁾= [ ˆA(qL, qR)1q]⁽ⁱ⁾, (17) for i = 1, 2, 3, where F ∈ R³is the usual definition of the fluxes for the Euler Eqs. (1), and1F⁽ⁱ⁾is the i th component of1F.

To accomplish the relation in (17), by taking a similar approach employed in [26] for real gases (cf. [23] also), we find it is sufficient to get the average states for variables such as ˆu, ˆH, ˆγ , ˆι, and set the matrix ˆA(qL, qR) = A(ˆu, ˆH, ˆγ , ˆι). The results are

ˆ

u =√ρ^LuL+ √ρRuR

√ρL+ √ρR

, Hˆ =√ρ^LHL+ √ρRHR

√ρL+ √ρR

,

(18) 1

γ − 1ˆ =√ρL

³ 1 γL−1

´+ √ρR

³ 1 γR−1

´

√ρL+ √ρR

, ˆι = √ρ^LιL+ √ρRιR

√ρL+ √ρR

.

Note that in the current derivation of the average states ˆγ and ˆι are chosen so that the expression

1p = ( ˆγ − 1)²

· 1

γ − 1ˆ 1ι − ˆι1 µ 1

γ − 1

¶¸

(11)

218 KEH-MING SHYUE

is satisfied approximately, and ˆp is computed by( ˆγ −1)ˆι (see [23, 24] for a related discussion to the case of real gases). As shown in Section 4, we find good results with the use of this set of average states defined in (18).

The solution of the linear problem (16) consists of five discontinuities propagating at constant speeds. 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 = X5 k=1

αˆkˆrk, (19)

where ˆrkis the kth eigenvector of ˆA with

ˆr1=





 1 ˆ u− ˆc Hˆ − ˆuˆc

0 0





, ˆr2=





 1 ˆ u

1 2uˆ²

0 0





, ˆr3=





 1 ˆ u+ ˆc Hˆ + ˆuˆc

0 0





, ˆr4=





 0 0 ˆp 1 0





, ˆr5=





 0 0 1 0 1





,

(20) and ˆA ˆrk= ˆλkˆrkwith ˆλkthe corresponding eigenvalue,

λˆ1= û − ˆc, λˆ2= û, λˆ3= û + ˆc, λˆ4= ˆλ5= û, (21) where ˆc=p

( ˆγ − 1)( ˆH − ˆu²/2) is the speed of sound. The scalar ˆαkgives the strength across the discontinuity that can be determined easily from (19). We find











αˆ2= γ − 1ˆ ˆc²

£( ˆH − ˆu²)1q⁽¹⁾+ ˆu1q⁽²⁾− 1q⁽³⁾+ ˆp1q⁽⁴⁾+ 1q⁽⁵⁾¤

= 1ρ − 1p ˆc² αˆ3= 1

2ˆc

£(ˆc − ˆu)1q⁽¹⁾+ 1q⁽²⁾− ˆc ˆα2

¤=1p + ˆρˆc1u 2ˆc² αˆ1= 1q⁽¹⁾− ˆα2− ˆα3= 1p − ˆρˆc1u

2ˆc² αˆ4= 1q⁽⁴⁾= 1

µ 1 γ − 1

¶

αˆ5= 1q⁽⁵⁾= 1 µγ p_∞

γ − 1

¶ ,

(22)

where ˆρ = √ρLρR.

Notice that in this Riemann solution there exist three discontinuities propagating at the same speed, ˆλ2= ˆλ4= ˆλ5= ˆu. For practical purposes, we may view these discontinuities as a single one with the operator ˆW2defined by combining all the jumps across the ˆλ2wave family, i.e., set ˆW2 = ˆα2ˆr2+ ˆα4ˆr4+ ˆα5ˆr5. Clearly, doing so removes the effect of the λˆ4 and ˆλ5wave families to the solution. With this notation, we also write ˆWk = ˆαkˆrkto represent the jump across the k= 1 and 3 waves. Wave propagation methods are based on using these propagating discontinuities to update the cells averages in the cells neighboring each interface.

We now show that in the case of an interface only problem the numerical solution obtained using the method would be free of oscillations, and in particular the pressure would remain

(12)

in equilibrium. Without loss of generality, we consider a single Riemann problem that at cell interface xjthe initial data are picked up so that there are no jumps in both the pressure and velocity, i.e.,1pⁿj = 0 and 1uⁿj = 0, but is otherwise for the other variables such as ρ, γ , and p∞. With this initial data, the Riemann solution consists of only a single contact discontinuity with the wave strength ˆαk jcomputed by (22) as

αˆ1 j= ˆα3 j = 0, αˆ2 j = 1ρⁿj, αˆ4 j = 1 µ 1

γ − 1

¶n

j

, αˆ5 j= 1 µγ p_∞

γ − 1

¶n

j

,

and the associated eigen-structure ˆrk j, ˆλk jcomputed by (20) and (21), respectively.

Now suppose that the wave moves to the right of the interface ˆuj > 0. To take account of the effect of this wave, according to (14), the cell average Uⁿ_j⁺¹should be evaluated by

ρⁿj⁺¹ = ρⁿj − 1t

1xuˆj1ρⁿj, (23a)

(ρu)ⁿj⁺¹ = (ρu)ⁿj− 1t

1xuˆj[ ˆu1ρⁿ]j, (23b)

(ρE)ⁿj⁺¹ = (ρE)ⁿj− 1t 1xuˆ_j

·1

2uˆ²1ρⁿ+ ˆp1 µ 1

γ − 1

¶n

+ 1 µγ p∞

γ − 1

¶n¸

j

, (23c) µ 1

γ − 1

¶n+1

j

= µ 1

γ − 1

¶n

j

− 1t 1xuˆj1

µ 1 γ − 1

¶n

j

, (23d)

µ γ p∞

γ − 1

¶n+1

j

= µγ p∞

γ − 1

¶n

j

− 1t 1xuˆj1

µγ p∞

γ − 1

¶n

j

. (23e)

It follows from (23a) and (23b) that we have the expected state of uⁿ_j⁺¹= uⁿjwhere by (18), ˆ

u_j = uⁿj. With this result, Eq. (23c) can be simplified to

(ρe)ⁿj⁺¹= (ρe)ⁿj− 1t 1xuˆ_j

· ˆp1

µ 1 γ − 1

¶n

+ 1 µγ p∞

γ − 1

¶n¸

j

,

or alternatively to µp+ γ p∞

γ − 1

¶n+1

j

=

µp+ γ p∞

γ − 1

¶n

j

− 1t 1xuˆ_j

· ˆp1

µ 1 γ − 1

¶n

+ 1 µγ p∞

γ − 1

¶n¸

j

,

when using the equation of state (2). Notice that in case (23e) is applied to the above equation, we have further simplification

µ p

γ − 1

¶n+1 j

=

µ p

γ − 1

¶n j

− 1t 1xuˆj

· ˆp1

µ 1 γ − 1

¶n¸

j

,

and find the pressure equilibrium of the computed solution pⁿ_j⁺¹ = pⁿj, when Eq. (23d) is satisfied along with the fulfillment of the condition ˆp_j= pⁿj (see (18) and the comment thereafter).

As to the behavior of the solutionsρⁿj⁺¹, γjⁿ⁺¹, and(p∞)ⁿj⁺¹, from (23a), (23d), (23e), it is easy to derive the modified equations for each of them individually, and show monotonicity

(13)

220 KEH-MING SHYUE

of the results obtained using the method (cf. [42]). See results shown in Section 4 for numerical verification of this statement also.

In practical applications, there may be some other waves which come into the j th cell and affect the cell average Uⁿ_j⁺¹as well. Suppose that we are taking Riemann data at cell interface x_j₊₁so that there is a 1-wave (shock or rarefaction) propagating to the left of the interface. In this method (14), we update the cell average Uⁿ_j⁺¹by

ρⁿj⁺¹ := ρⁿj⁺¹− 1t

1x(ˆu − ˆc)j+1

·1pⁿ− ˆρˆc1uⁿ 2ˆc²

¸

j+1, (ρu)ⁿj⁺¹ := (ρu)ⁿj⁺¹− 1t

1x(ˆu − ˆc)j+1

· (ˆu − ˆc)

µ1pⁿ− ˆρˆc1uⁿ 2ˆc²

¶¸

j+1, (ρE)ⁿj⁺¹ := (ρE)ⁿj⁺¹− 1t

1x(ˆu − ˆc)j+1

·

( ˆH − ˆuˆc)

µ1pⁿ− ˆρˆc1uⁿ 2ˆc²

¶¸

j+1. We note that due to a fundamental property of the Roe solver, i.e., the relation in (17), this is a conservative update of the numerical solutionsρⁿj⁺¹, (ρu)ⁿj⁺¹, and(ρE)ⁿj⁺¹. We get a new pressure pⁿ_j⁺¹by using the equation of state (2),

pⁿ_j⁺¹=¡

γjⁿ⁺¹− 1¢µ ρE −1

2ρu²

¶n+1

j

− (γ p∞)ⁿj⁺¹.

One advantage of using the wave propagation form is that we are able to handle each wave in turn, and there is no need to compute fluxes and make distinction between the waves (cf. [67]). Extension of the method to higher order accuracy, and in particular to a high-resolution version of the wave propagation scheme, follows easily by incorporating limited slopes and by constructing piecewise linear profiles to the method; see [43, 44]

for the detail. It is not difficult to show that for the interface only problem we again have the required pressure equilibrium that is independent of the limiter being employed to the high-resolution method. Moreover, we obtain a better resolution of the result as compared to the first order result; see Section 4 for numerical examples.

3.2. Method with volume-fraction model. We now consider the use of the Riemann solutions associated with the volume-fraction model (13) to the wave propagation method (14). For simplicity in description, we are concerned with a two-component flow problem that the quasi-linear system of equations in (15) are defined with

q=





 ρ ρu ρE Y⁽¹⁾ Y⁽²⁾





, A(q) =







0 1 0 0 0

³γ −3 2

´

u² (3 − γ )u γ − 1 χ⁽¹⁾ χ⁽²⁾

³_{γ −1}

2

´

u³− uH H − (γ − 1)u² γ u χ⁽¹⁾u χ⁽²⁾u

0 0 0 u 0

0 0 0 0 u





 ,

whereχ⁽ⁱ⁾= (1 − γ )(p + γ⁽ⁱ⁾p⁽ⁱ⁾_∞)/(γ⁽ⁱ⁾− 1). We use the Roe solver as usual in that we solve the linear problem (16) with the Roe-type matrix ˆA(qL, qR) = A(ˆu, ˆH, ˆι, ˆY⁽¹⁾, ˆY⁽²⁾) depending on the initial data qL and qR. Here as in the case of the γ -based model we

(14)

compute the average states ˆu, ˆH, ˆι by (18), and set the average state ˆY⁽ⁱ⁾by the standard

“Roe-averaging” approach also with

Yˆ⁽ⁱ⁾= √ρ^LY_L⁽ⁱ⁾+ √ρRY_R⁽ⁱ⁾

√ρL+ √ρR

for i = 1, 2.

When defining ˆA in this way, it is easy to check the satisfaction of the fundamental relation in (17).

The result of this linearized Riemann problem takes a rather similar wave structure as to the form shown in (19)–(22) of theγ -based model. Without writing the full solution here, we mention that the only difference between the representation of these two solutions is to the jump across the ˆλ4and ˆλ5families. Here we use the following expressions instead,

αˆ4= 1Y⁽¹⁾, αˆ5 = 1Y⁽²⁾, ˆr4 =





 0 0

ˆ p+ γ⁽¹⁾p⁽¹⁾_∞

γ⁽¹⁾− 1

1 0





, ˆr5=





 0 0

ˆp+ γ⁽²⁾p⁽²⁾_∞ γ⁽²⁾− 1

0 1





.

As in the case performed for theγ -based model, we now analyze the pressure obtained using the method (14) with the volume-fraction model for the interface only problem con- sidered there. In this case, it is enough to begin looking at the update of the internal energy in the form,

µp+ γ p_∞ γ − 1

¶n+1 j

=

µp+ γ p_∞ γ − 1

¶n j

− 1t 1xuˆj

" ₂ X

i=1

ˆp+ γ⁽ⁱ⁾p⁽ⁱ⁾_∞ γ⁽ⁱ⁾− 1 1¡

Y⁽ⁱ⁾¢n

#

j

.

Note the above equation is just a simplified version of the update of the total energy by taking into account of the fact that uⁿ_j⁺¹= uⁿj. When substituting the volume-fraction relations in (6a) and (6c) forγ and p_∞, we may write the equation as

Ã ₂ X

i=1

Y⁽ⁱ⁾p

γ⁽ⁱ⁾− 1+Y⁽ⁱ⁾γ⁽ⁱ⁾p_∞⁽ⁱ⁾ γ⁽ⁱ⁾− 1

!n+1

j

= Ã ₂

X

i=1

Y⁽ⁱ⁾p

γ⁽ⁱ⁾− 1 +Y⁽ⁱ⁾γ⁽ⁱ⁾p⁽ⁱ⁾_∞ γ⁽ⁱ⁾− 1

!n

j

−1t 1xuˆj

" ₂ X

i=1

ˆp+ γ⁽ⁱ⁾p⁽ⁱ⁾_∞ γ⁽ⁱ⁾− 1 1¡

Y⁽ⁱ⁾¢n

#

j

.

After a simple manipulation, we find pⁿ_j⁺¹= pⁿj, only if ˆp_j = pⁿjand the establishment of the difference equations for the volume fraction Y⁽ⁱ⁾,

¡Y⁽ⁱ⁾¢n+1

j =¡

Y⁽ⁱ⁾¢n j− 1t

1xuˆ_j1¡ Y⁽ⁱ⁾¢n

j for i = 1, 2.

It is easy to see that the latter two requirements are guaranteed by the method.

Note for this two-component flow problem we have described the method uses two different volume-fraction functions Y⁽¹⁾and Y⁽²⁾. In practice, due to the fact that Y⁽¹⁾+ Y⁽²⁾= 1, we may simply use a volume-fraction function, say Y⁽¹⁾, to the model, and set the

(15)

222 KEH-MING SHYUE

value of the other one, Y⁽²⁾= 1 − Y⁽¹⁾, for example. When modeling in this way, it gives a more robust method to the computation of two-component flow problems than the use of theγ -based model to the method.

4. NUMERICAL RESULTS IN ONE DIMENSION

We now present results to validate our multicomponent algorithm with the Roe solver described in Section 3 in one dimension.

EXAMPLE4.1. We consider an interface only problem that the solution of a Riemann problem consists of a single contact discontinuity in gas dynamics. We take two sets of data for numerical experiments. In the first set, we have a polytropic gas and use two constant states as





 ρ u p γ p_∞







L

=





 1 1 1 1.4

0





 and





 ρ u p γ p_∞







R

=





 0.125

1 1 1.2

0





, (24)

while in the second set we use the same data as in the first set with the exception that on the state R a stiffened gas withγR = 4 and (p∞)R = 1 is employed instead. Here L is the state used for x∈ [0, 0.2) and R is the state used for x ∈ [0.2, 1].

Results for the first set of data are shown in Fig. 1 where we have employed theγ -based model together with both the first order and high-resolution wave propagation methods to the computations. For comparison, we also include results obtained using three different effective equations to the simulations. They are the conservation equations with either

∂

∂t(ργ ) + ∂

∂x(ργ u) = 0, (25a)

or

∂

∂t µ ρ

γ − 1

¶ + ∂

∂x µ ρu

γ − 1

¶

= 0, (25b)

and the primitive equation

∂γ

∂t + u∂γ

∂x = 0. (25c)

Note in the cases of (25a) and (25b), we have two model systems that are in the full conservation form; see [34] for a similar consideration.

From the figure, we clearly observe pressure fluctuations in the solutions when employing any of the equations in (25a)–(25c), but not theγ -based model where (9) is considered, to the method. By following a similar analysis conducted in Section 3 for an interface only problem, we may explain the observed error behavior in pressure as being the failure to approximate the energy Eq. (7b) in a consistent manner when those equations are in use (cf. [34]). Note other variables in the solution are affected by this error also as time evolves.

Numerical evidences suggest however that these errors decrease as the mesh is refined, and

(16)

FIG. 1. Comparison plots of four different multicomponent models for the Euler equations with data (24) at time t = 0.12. (a) Results using the first order wave propagation method. (b) Results using the high-resolution wave propagation method with the “minmod” slope limiter. In each figure the solid line is the exact solution and the points show the computed solution with 100 mesh points.

(17)

224 KEH-MING SHYUE

FIG. 2. High-resolution results for an interface only problem using the same initial data as in Fig. 1 with the exception thatγR= 4 and (p∞)R= 1 are used in the current computation. The solid line is the exact solution and the points show the computed solution with 100 mesh points.

the rate of convergence of the error in the 1-norm is about the order of accuracy of the method that is employed to the computation.

Figure 2 shows results of a run using the second set of data where only the high-resolution solutions with the γ -based model are presented. Notice that the pressure and also the velocity remain at constant states for this stiffened gas simulation. Because the errors become too erroneous in most cases, we do not show results using the type of models given in (25a)–(25c) for this test.

Note in the above computation, we use 100 mesh points and plot the results at time t= 0.12. For the purposes in illustration of the basic solution structure obtained from using the high-resolution version of the method, we only present results using the sim- pler “minmod” slope limiter (see [42, 70]) for the runs. Of course, other more sophisticated limiters, such as “superbee” for example (cf. [2, 70]), can be employed to the methods for computations also. As far as the global structure of the solution is concerned, we observe quite similar behavior of the solutions when different limiters are in use to the method. We further remark that the Courant number (see [42]) in choosing the time step that maintains stability of the method is 0.9 in the tests. The nonreflecting boundaries are used on the left and right of the computational domain. Without further notice, we use the same limiter and the Courant number for all other experiments performed in this paper.

EXAMPLE4.2. We next are concerned with a two-phase gas-liquid Riemann problem.

On the left when x∈ [0, 0.5), we have the gas phase with data (ρ, u, p, γ, p_∞)L = (1.241, 0, 2.753, 1.4, 0), and on the right when x∈ [0.5, 1], we have the liquid phase with data

(ρ, u, p, γ, p∞)R= (0.991, 0, 3.059 × 10⁻⁴, 5.5, 1.505).

We note that the above variables have been nondimensionalized as in the work done by Cocchi et al. [12] and Cooke and Chen [17] to simulate underwater explosions in a spheri- cally symmetric geometry. We run the problem in a shock tube with 100 mesh points, and show the high-resolution results in Fig. 3 at time t = 0.1 using the γ -based model. From the figure, we again see the correct behavior of the computed contact discontinuity, and

(18)

FIG. 3. High-resolution results for a two-phase gas-liquid Riemann problem at time t= 0.1. The solid line is the exact solution and the points show the computed solution with 100 mesh points.

also the rarefaction and shock waves as in comparison with the exact solution (the solid line shown in the plot). The results of this problem with a source term that accounts for the simplification of a cylindrically symmetric flow will be presented in Section 5 together with two-dimensional results.

EXAMPLE 4.3. Finally, we consider a shock-contact interaction problem studied by Abgrall [1] and Karni [35] that verifies convergence of the computed solutions to the correct weak ones in a multicomponent case. The initial condition we use consists of a stationary interface at x = 0.5 separating two fluids of different equation of states, and a leftward going Mach 1.95 shock wave at x = 0.6 traveling from the right to left. The gas on the left of the interface is a polytropic gas with

(ρ, u, p, γ, p∞)L = (1, 0, 1, 1.4, 0),

and the gas on the right of the interface (i.e., on the middle and the preshock state), is a stiffened gas with

(ρ, u, p, γ, p_∞)M = (5, 0, 1, 4, 1).

The state behind the shock is

(ρ, u, p, γ, p∞)R = (7.093, −0.7288, 10, 4, 1).

The exact solution for this problem in the x-t plane up to time t= 0.2 is illustrated in Fig. 4 where the density contours are presented.

A snap shot of the computed total internal energy, velocity, and pressure are shown in Fig. 5 at time t= 0.2 where we solve the problem using the high-resolution wave propagation method with 200 mesh points. We can easily see that the shock wave and contact discontinuity are very well located, and the rarefaction wave moves at the correct speed with the correct shape. A two-dimensional version of the problem will be considered in Section 6.

We note that in this section we have only present numerical solutions obtained using theγ -based model to the method. It is interesting to mention that we find little difference between the results as compared to the ones obtained using the volume-fraction model to the method for simulations. Because of this, we omit the presentation of the volume-fraction based numerical results here.