The coupling of homogeneous models for two-phase...

The coupling of homogeneous models for two-phaseflows

Annalisa Ambroso1, Christophe Chalons2, Frederic Coquel3, Edwige Godlewski3,Frederic Lagoutiere2, Pierre-Arnaud Raviart3, Nicolas Seguin3

1 DEN/DM2S/SFME/LETR CEA-Saclay, F-91191 Gif-sur-Yvette, France,2 Universite Paris 7-Denis Diderot et UMR 7598 Laboratoire Jacques-Louis Lions, Paris,

F-75005 France,3 Universite Pierre et Marie Curie-Paris6, UMR 7598 LJLL, Paris, F-75005 France;

CNRS, UMR 7598 LJLL, Paris, F-75005 France

[email protected], [email protected], [email protected],

[email protected], [email protected], [email protected],

[email protected]

Abstract

We consider the numerical coupling at a fixed spatial interface of two

homogeneous models used for describing non isothermal compressible two

phase flows. More precisely, we concentrate on the numerical coupling

of the homogeneous equilibrium model and the homogeneous relaxation

model in the context of finite volume methods. Three methods of coupling

are presented. They are based on one of the following requirements:

continuity of the conservative variable through the coupling interface,

continuity of the primitive variable and global conservation of mass,

momentum and energy. At the end, several numerical experiments are

presented in order to illustrate the ability of each method to provide

results in agreement with their principle of construction.

1 Introduction

The coupling of fluid flow models is becoming a key topic in industrial code develop-ment. From an engineering point of view, different models are used to treat differentsub-domains of a complex system where a flow takes place. Therefore it is naturalto want to have a global simulation for the system by putting these domains, andthus these models, side to side. This rises the delicate question of the continuity ofthe description of the flow. As an example, think of the simulation of a combustionengine, where the natural models for the fuel pipes, the injector and the combustionchamber are clearly different, or, again, think of the treatment, in nuclear energyindustry, of the coolant circuits which are formed by different components, each onewith its associated specific model for the coolant flow. Our work is motivated bythis last application and the coolant under consideration is a two-phase fluid.

The coupling of homogeneous models for two-phase flows

Two-phase flows can be described by means of different models: mixture, drift(homogeneous or not), two-fluid or even multi-field models are currently used inindustrial thermo-hydraulic codes. We consider here the problem of the interfacialcoupling of two models, i.e. we imagine that a two-phase flow is described by meansof a model M1 at the left of a fictitious interface I and by another model M2 atthe right of I. Our aim is to numerically describe the whole flow dealing with thepotential discontinuities at I, when the model jumps. It is clear that this analysisstrongly depends on the two models M1 and M2 that describe the flow on each sideof the discontinuity.

In this paper we focus on the coupling of two homogeneous models. More pre-cisely, we consider as M1 a homogeneous equilibrium model frequently referred toas HEM and, as model M2 a homogeneous relaxation model HRM. Homogeneousmodels describe a two-component flow as the flow of a single fluid: the compound.A description of these models can be found in [10], [22]. When the two componentsare assumed to follow a perfect gas law, the full thermodynamic of the compoundcan be described by means of analytic formulas, as it is clearly presented in [30] andin [13], [15], [16], [17]. We will focus on this case and, for the sake of completeness,we present these models in Section 2. The HEM model is detailed in Section 2.2.1while Section 2.2.2 details the HRM one.

The general framework of the coupling problem is as follows. Let D be anopen subset of R

n (n ∈ N, n ≥ 1) divided in two separated open sub-domains DE

and DR and an interface I: D = DE ∪ DR ∪ I. We assume that the flow evolvesaccording to HEM in DE and HRM in DR. We note uE and uR the vectors of thecorresponding conservative unknowns. In the following, we restrict the study to theone-dimensional case where DE = R

−, DR = R+ and I = 0. Therefore, the flow

is described by the following problem

HEM on R+ ×DE

∂tuE + ∂xf

E(uE) = SE(uE),

∣∣∣∣∣∣∣∣∣∣∣∣∣∣

HRM on R+ ×DR

∂tuR + ∂xf

R(uR) = SR(uR),

(1)

where fE, fR and SE,SR respectively denote the conservative flux and the sourceterms of the models, that will be given in Section 2. We underline the lack ofinformation at the interface I. Our only constraint is given by the context of theproblem. We know that the physical entity we are describing is the same: the flowof the same fluid. Its properties are characterized by uE on one side of the interfaceand by uR on the other side, and we wish to link these quantities at the crossing ofI taking into account the continuity of the flow.

We acknowledge the artificial character of this problem: the real physics is thesame, but, when different codes are used to simulate it, the global description can bediscontinuous. This can be seen as an artificial error introduced by the simulationprocess. Our analysis wants to deal with this error and its effects on the solution ofthe problem.

International Journal on Finite Volumes 2


First of all, the problem must be analyzed from the modelling point of view.To make (1) complete, we need to add an interface model, that is to say, we mustdetail the properties of the flow that should hold at the crossing of the interface. Wethink, for instance, of the continuity of some variables or the conservation of somequantities. In a second place, the global problem should be analyzed mathemati-cally. Compatibility of the constraints imposed at the interface within one anotherand with the models should be checked. Moreover, we point out that, due to thehyperbolic character of the models under consideration (see Section 2), it is verylikely that a given coupling problem with a prescribed model for the interface can-not be solved for all initial data. Finally, we must not forget that we are interestedin the numerical simulation of the flow. This has two main implications. On theone hand, we wish that the constraints given by the chosen interface model leadto simple numerical relations. On the other hand, we must keep in mind that thestrategy of coupling we choose can be effective only if its numerical treatment isappropriate.

The problem of coupling nonlinear conservation laws was first addressed in [25]and in [26]. In [4] and in [5] a study of the case where the models to be coupled aregiven by two Euler systems for gas dynamics with two different equations of state ispresented. In this case, one of the numerical techniques proposed is very similar tothe ghost fluid method ([32], [2]). The coupling between a one-dimensional model ofgas dynamics and a 2-dimensional one are treated in [28]. References [27] and [11]study the case of a change of porosity of the medium where the flow takes place,while [33] and [35] deal with the coupling of acoustic wave equations with differentsound speed. For the case under study in the present article, i.e. the coupling ofa Homogeneous Relaxation Model with a Homogeneous Equilibrium Model, someresults were presented in [6] and a study of an industrial case of interest can befound in [29].

In this paper, we propose several ways for numerically solving the coupling prob-lem (1), depending on the treatment we impose at the interface I. The HRM andHEM models are presented in Section 2. The interface models we consider are listedin Section 3. In Section 4 we give the basis for the numerical treatment of theproblem. In particular, we present the two main strategies that can be chosen: theflux coupling (Section 4.1.1) and the intermediate state coupling (Section 4.1.2).Section 4 ends with a description of the numerical schemes we consider.Numerical tests to illustrate the different choices are presented in Section 5, followedby a conclusion (Section 6).

2 Governing equations and physical modelling

In this study, we focus on the one-dimensional case (n = 1) and we set D = R,DE = R

−,∗ and DR = R+,∗. The coupling interface is then located at x = 0:

I = x = 0. Let us now precise the general modelling assumptions on the two-phase flow.



2.1 General modelling assumptions

We follow [13], [15], [16], [17], [30] and consider that each phase is a fluid with itsown thermodynamic properties. Let us begin by introducing some notations for thetwo phases. In the following, ρα, εα, pα(ρα, εα) and Tα(ρα, εα) respectively denotethe density, the internal energy, the pressure and the temperature of the phaseα = 1, 2. Then, the entropy ηα(ρα, εα) is defined up to an additive constant viaGibbs relations, that is

∂ηα

∂εα(ρα, εα) =

1

Tα(ρα, εα),

∂ηα

∂ρα(ρα, εα) = −

pα(ρα, εα)

ρ2αTα(ρα, εα)

. (2)

In order to make the localization of phase changes possible, we introduce an orderparameter z ∈ [0, 1] that equals 1 in phase 1 and 0 in phase 2.

The matter is now to put these two fluids in relation in a way to obtain athermodynamically coherent description of the two-phase flow. With this in mind,we adopt the following definitions for the mixture density ρ, internal energy ε andpressure p:

ρ = zρ1 + (1 − z)ρ2,ρε = zρ1ε1 + (1 − z)ρ2ε2,p = zp1(ρ1, ε1) + (1 − z)p2(ρ2, ε2).

(3)

Then, we make the assumption that the two fluids have the same temperature.Therefore, we impose the following closure relation which defines the mixture tem-perature T

T = T1(ρ1, ε1) = T2(ρ2, ε2). (4)

We complete this section by giving explicit formulas for pα and Tα, α = 1, 2. Weconsider the case of perfect gas equations of state given by

pα(ρα, εα) = (γα − 1)ραεα,

Tα(ρα, εα) =εα

Cv,α, α = 1, 2, (5)

with different adiabatic coefficients γ1 > 1 and γ2 > 1. We choose γ2 < γ1 withoutloss of generality. For simplicity reasons, specific heats Cv,1 and Cv,2 are assumedto be equal and we set

Cv = Cv,1 = Cv,2 > 0. (6)

This hypothesis is widely used in the litterature ([13], [15], [16], [17], [30]) and allowsto ease the computations. For the general case with different specific heat coefficientsfor the two fluids, we refer to [23] and [3].As an immediate consequence of (3)-(4)-(6), we thus have

ε = ε1 = ε2 = CvT. (7)

Integrating Gibbs relations (2), we then define entropies ηα by

ηα(ρα, εα) = Cv,α lnpα(ρα, εα)

ργαα

, α = 1, 2. (8)



2.2 Governing equations

In this section, we describe the basic assumptions that lead to the homogeneousequilibrium model HEM and to the homogeneous relaxation model HRM and wegive the two sets of equations associated with them. We begin by HEM.

2.2.1 The homogeneous equilibrium model HEM

The homogeneous equilibrium model is obtained under specific assumptions of me-chanical and thermodynamic equilibrium for the two phase flow. This means thatwhen both phases are present, their pressures pα(ρα, εα) and their free enthalpiesgα(ρα, εα) = εα + pα(ρα, εα)/ρα − ηα(ρα, εα)Tα(ρα, εα) are equal:

p1(ρ1, ε1) = p2(ρ2, ε2),g1(ρ1, ε1) = g2(ρ2, ε2),

(9)

For the perfect gas mixture described above, this is equivalent to

(γ1 − 1)ρ1 = (γ2 − 1)ρ2,

(γ1 − 1) − ln(γ1 − 1)

ρ(γ1−1)1

= (γ2 − 1) − ln(γ2 − 1)

ρ(γ2−1)2

, (10)

thanks to (5)-(7)-(8). By easy calculations, we get

Proposition 2.1 System (10) admits a unique solution (ρ?1, ρ

?2) given by

ρ?1 =

1

e(γ2 − 1

γ1 − 1)

γ2

γ2 − γ1 ,

ρ?2 =

1

e(γ2 − 1

γ1 − 1)

γ1

γ2 − γ1 .

(11)

Remark 1 Under assumption γ2 < γ1, we have ρ?1 < ρ?

2.

At this stage, it is important to notice that considering perfect gas equations ofstate with the same specific heats leads to explicit and simple formulas for equilib-rium densities ρ?

1 and ρ?2. This is actually the main motivation for such a choice.

As an immediate consequence of (6), note also that ρ?1 and ρ?

2 do not depend on themixture temperature T .

Remark 2 Entropies ηα being defined up to an additive constant, one could have set

ηα(ρα, εα) = Cv,α lnpα(ρα, εα)

ργαα

− γαCv,α ln(γα − 1)Cv,α, α = 1, 2, (12)

instead of (8). In this case, ρ?1 and ρ?

2 are given by

ρ?1 = exp(−1 −

γ2 ln(γ2 − 1)Cv − γ1 ln(γ1 − 1)Cv

γ2 − γ1) × (

γ2 − 1

γ1 − 1)

γ2

γ2 − γ1 ,

ρ?2 = exp(−1 −

γ2 ln(γ2 − 1)Cv − γ1 ln(γ1 − 1)Cv

γ2 − γ1) × (

γ2 − 1

γ1 − 1)

γ1

γ2 − γ1 ,

(13)



or, in an equivalent way,

ρ?1 = exp

(γ1 + (γ1 − 1) ln(γ1 − 1)Cv − γ2 − (γ2 − 1) ln(γ2 − 1)Cv

γ2 − γ1

)

×

(γ2 − 1

γ1 − 1

) γ2 − 1

γ2 − γ1 ,

ρ?2 = exp

(γ1 + (γ1 − 1) ln(γ1 − 1)Cv − γ2 − (γ2 − 1) ln(γ2 − 1)Cv

γ2 − γ1

)

×

(γ2 − 1

γ1 − 1

) γ1 − 1

γ2 − γ1 ,

(14)

which corresponds to the values proposed in [13].

This result implies that in the HEM framework, the densities ρα of each phaseare given fixed values when both phases are present, i.e. when order parameterz ∈ (0, 1). The mixture density ρ = ρ?

1z + ρ?2(1 − z) thus lies in the interval

[ρ?1, ρ

?2] and depends on the value of z. Concerning the pressure law p, we have

p(ρ, ε) = p1(ρ?1, ε1) = p2(ρ

?2, ε2) for all ρ ∈ [ρ?

1, ρ?2] which corresponds to the first

assumption in (9). Outside this interval, it can be naturally extended since onlyphase 1 (respectively phase 2) is present when 0 < ρ ≤ ρ?

1 (respectively ρ ≥ ρ?2). We

get

p(ρ, ε) =

(γ1 − 1)ρε if ρ ≤ ρ?1,

(γ1 − 1)ρ?1ε = (γ2 − 1)Cvρ

?2T if ρ?

1 < ρ < ρ?2,

(γ2 − 1)ρε if ρ ≥ ρ?2.

(15)

In the following, this pressure law will be noted pE(ρ, ε).Then, the homogeneous equilibrium model describing the one-dimensional flow

under consideration is given by

∂tρ + ∂x(ρu) = 0,∂t(ρu) + ∂x(ρu2 + p) = 0,∂t(ρE) + ∂x(ρE + p)u = 0.

(16)

This corresponds to uE = (ρ, ρu, ρE), fE(uE) = (ρu, ρu2+p, ρEu+pu) and SE(uE) =(0, 0, 0) with the notations introduced in (1). The first equation expresses the conser-vation of mass, while the second and the third respectively govern the conservationof momentum ρu and total energy ρE. Note that the two fluids share the same ve-locity u which is inherent in the homogeneous modelling. The pressure p = pE(ρ, ε)is given by (15) while internal energy ε is linked to the vector uE by

ρE = ρε +1

2ρu2. (17)

Defining the following natural phase space for HEM:

ΩE = uE := (ρ, ρu, ρE) ∈ R3/ρ > 0, ε = CvT > 0,

we are in position to state:



Lemma 2.2 The first order convective system HEM is strictly hyperbolic over ΩE,with the following eigenvalues:

λ1(uE) = u − cE < λ2(u

E) = u < λ3(uE) = u + cE,

where the sound speed cE is such that

(cE(uE)

)2

=

γ1(γ1 − 1)ε if 0 < ρ ≤ ρ?1,

(γ1 − 1)2ρ?12

ρ2ε = (γ2 − 1)2

ρ?22

ρ2ε if ρ?

1 < ρ < ρ?2,

γ2(γ2 − 1)ε if ρ ≥ ρ?2.

(18)

Moreover, the second field is linearly degenerate.

The proof of this lemma follows from standard calculations and is left to the reader.See also [30], [13], [15], [16], [17].

Remark 3 It is worth noting that the sound speed uE → cE(uE) is discontinuouswhen the density ρ is equal to ρ?

1 or ρ?2. As a result, the first and the third fields,

corresponding respectively to the eigenvalues u−cE and u+cE, are neither genuinelynonlinear nor linearly degenerate. This property is known to make more complicatedthe resolution of the Riemann problem since admissible nonclassical solutions vio-lating the standard selection criterion naturally arise. Such considerations are ofcourse out of the scope of this paper. Besides, existence and uniqueness may berecovered when imposing the validity of the Liu criterion for instance. We refer forinstance the reader to [30], [31] and the references therein.

Let us now address the homogeneous relaxation model.

2.2.2 The homogeneous relaxation model HRM

The homogeneous relaxation model considers that the two-phase flow no longerevolves instantaneously at thermodynamic equilibrium, but only at mechanical equi-librium. Modelling assumptions (9) are replaced by

p1(ρ1, ε1) = p2(ρ2, ε2), (19)

that is, in the case under study,

(γ1 − 1)ρ1 = (γ2 − 1)ρ2 (20)

due to (5) and (7). In other words, and contrary to HEM, densities ρ1 and ρ2 are notrestricted any longer to take the saturation values ρ?

1 and ρ?2, but are simply linked

by the algebraic relation (20). Actually, HRM accounts for mass transfers betweenthe two fluids assuming that the thermodynamic equilibrium g1(ρ1, ε1) = g2(ρ2, ε2)is not instantaneously achieved, but it is reached at speed λ0 > 0. More precisely,the system reads

∂t(ρ1z) + ∂x(ρ1zu) = λ0

(ρ?1z

?(ρ) − ρ1z),

∂tρ + ∂x(ρu) = 0,∂t(ρu) + ∂x(ρu2 + p) = 0,∂t(ρE) + ∂x(ρE + p)u = 0,

(21)



that is, with the notations introduced in (1), uR = (ρ1z, ρ, ρu, ρE), fR(uR) =(ρ1zu, ρu, ρu2 + p, ρEu + pu) and SR(uR) = (λ0

(ρ?1z

?(ρ) − ρ1z), 0, 0, 0). In order

to close system (21), let us recall that ρ = ρ1z + ρ2(1 − z), therefore

z =ρ − ρ2

ρ1 − ρ2,

while z?(ρ) corresponds to the thermodynamic equilibrium value given by ρ =ρ?1z

?(ρ) + ρ?2(1− z?(ρ)) when ρ ∈ [ρ?

1, ρ?2]. Otherwise, we naturally set ρ?

1z?(ρ) = ρ if

ρ ≤ ρ?1 and ρ?

1z?(ρ) = 0 if ρ ≥ ρ?

2, so that

z?(ρ) =

ρ

ρ?1

if 0 < ρ ≤ ρ?1,

ρ − ρ?2

ρ?1 − ρ?

2

if ρ?1 ≤ ρ ≤ ρ?

2,

0 if ρ ≥ ρ?2.

(22)

Note that ρ2 > ρ1 since we have assumed γ1 > γ2. Pressure p simply follows fromidentity p = zp1(ρ1, ε) + (1 − z)p2(ρ2, ε), that is

p = p(ρ1z, ρ, ε) =((γ1 − 1)ρ1z + (γ2 − 1)(ρ − ρ1z)

)ε, (23)

and (17) remains valid. Hereafter, this pressure law will be noted pR(ρ1z, ρ, ε).

Remark 4 We note that in the limit λ0 → +∞, usually called equilibrium, HRMconverges at least formally toward HEM. Indeed, the first equation in (21) leads inthis asymptotic regime to the relation ρ1z = ρ?

1z?(ρ) and we get the pressure law

(15) by means of (23) and (22). In other words, pR,eq := pR(ρ?1z

?(ρ), ρ, ε)

equalspE := pE(ρ, ε).

To conclude, a natural phase space for HRM is

ΩR = uR := (ρ1z, ρ, ρu, ρE) ∈ R4/ρ > 0, 0 ≤ ρ1z ≤ ρ, ε = CvT > 0,

and the following statement holds true.

Lemma 2.3 The first order underlying system of HRM is hyperbolic over ΩR, withthe following eigenvalues:

λ1(uR) = u − cR < λ2(u

R) = λ3(uR) = u < λ4(u

R) = u + cR,

where the sound speed cR is such that

(cR(uR)

)2

=A(ρ1z, ρ)

ρ

(1+

A(ρ1z, ρ)

ρ

)ε, A(ρ1z, ρ) = (γ1−1)ρ1z+(γ2−1)(ρ−ρ1z).

(24)Moreover, the first and fourth fields are genuinely nonlinear and the second and thethird ones are linearly degenerate.



Remark 5 Let uE be a vector of ΩE and uR,eq =(ρ?1z

?(ρ),uE)

the associated vec-tor of ΩR taken at equilibrium. Then, it can be shown by easy calculations thatcR(uR,eq) ≥ cE(uE) which ensures that the first (respectively last) eigenvalue of themodel HRM is less or equal (respectively greater or equal) than the first (respectivelylast) eigenvalue of the model HEM. This property can be related to the so-called sub-characteristic condition (or Whitham condition) that guarantees the stability of therelaxation process in relaxation systems. See for instance [9], [18] and the literatureon this subject.

3 The coupling problem

The coupling problem we study is as follows. The global problem is 1D (see [28] forcoupling problems in higher dimensions). The interface coupling is I = x = 0 andseparates the two open sub-domains DE = R

−,∗ and DR = R+,∗. The flow is thus

governed by the following set of PDE:

∂tuE + ∂xf

E(uE) = SE(uE), for t > 0 and x ∈ DE,∂tu

R + ∂xfR(uR) = SR(uR), for t > 0 and x ∈ DR,

detailed respectively in (16) and (21), complemented by the initial data

uE(x, 0) = uE0 (x), for x ∈ DE,

uR(x, 0) = uR0 (x), for x ∈ DR,

where uE0 ∈ ΩE and uR

0 ∈ ΩR are given.It remains to define the behavior of the flow through I. Since the two models are

different, there is no obvious way to couple the systems HEM and HRM. However,these models are somehow compatible (see Rem. 4) and govern the evolution ofthe same fluid. Therefore, the coupling should be performed in such a way thatthe solution obeys some physical requirements. Here, we provide three numericalcoupling methods, in order to ensure respectively:

• the global conservation of ρ, ρu and ρE,

• the continuity of ρ, ρu and ρE (and the conservation of ρ) through I,

• the continuity of ρ, ρu and p (and the conservation of ρ and ρu) through I.

(For rigorous definitions of coupling, we refer to [25, 26].) Let us emphasize that thetwo latter coupling must be understood in a weak sense, since we are dealing withhyperbolic systems and discontinuous solutions. Indeed, if the characteristic speedsare incompatible, the variables expected to be constant through I can jump. As anexample, think about the transport equation ∂tv + a∂xv = 0, with a > 0 for x < 0and a < 0 for x > 0. In this case, the solution will be discontinuous [25].

In the following, a numerical method of coupling is proposed for each interfacecoupling. Afterward, the ability of each numerical coupling method to provide ap-proximate solutions in agreement with the corresponding coupling condition (that iseither the global conservation, the continuity of the conservative variable or the con-tinuity of the primitive variable) is illustrated by numerous numerical experiments.



4 The numerical coupling

In this section, we show how to couple the two homogeneous models from a numericalpoint of view. We aim at presenting several strategies, depending on informationsthat should be transferred through the interface and/or expected properties of thesolution, in terms of conservation in particular.

Let be given a constant time step ∆t and a constant space step ∆x and let usset ν = ∆t/∆x. Introducing the cell interfaces xj = j∆x for j ∈ Z and tn = n∆t forn ∈ N, we classically seek at each time tn a piecewise constant approximate solutionx → uν(x, tn) of the solution u of the coupling problem (1):

uν(x, tn) = unj+1/2 for all x ∈ Cj+1/2 = [xj;xj+1), j ∈ Z, n ∈ N.

Note that for j = 0, xj = x0 coincides with the coupling interface and, for alln ∈ N, un

j+1/2 = (ρ, ρu, ρE)nj+1/2 has three components for j < 0 and un

j+1/2 =

(ρ1z, ρ, ρu, ρE)nj+1/2 has four components for j ≥ 0.

In order to advance some given sequence (unj+1/2)j∈Z at time tn to the next time

level tn+1, our approach proposes two steps based on a splitting strategy betweenthe convective parts and the source terms of HEM and HRM. More precisely:

Step 1 (tn → tn+1−)In the first step, we focus on the convective part of the coupling problem. It consistsin solving

∂tuE + ∂xf

E(uE) = 0, x < 0,∂tu

R + ∂xfR(uR) = 0, x > 0,

(25)

for t ∈ (0,∆t] with initial data uν(., tn) and some coupling condition which we

describe below. We denote uν(., tn+1−) the corresponding solution at time t = ∆t.

Step 2 (tn+1− → tn+1)Then, uν(., t

n+1−) naturally serves as initial data for solving the source terms:

∂tuE = SE(uE), x < 0,

∂tuR = SR(uR), x > 0,

(26)

again for t ∈ (0,∆t]. This step will provide us with an updated solution uν(., tn+1)

and this completes the numerical procedure.

4.1 The convective part

Let us first address the first step devoted to the convective part. In this context, weassume two schemes to be given, by mean of two 2-point numerical flux functionsgE : R

3 × R3 7→ R

3 and gR : R4 × R

4 7→ R4 respectively consistent with the flux

functions fE and fR in the sense of finite volume methods. The former are used tonumerically solve the first-order underlying systems of HEM and HRM on each sideof the interface, that is

un+1−j−1/2 = un

j−1/2 − ν(gEj − gE

j−1), j ≤ 0,

un+1−j+1/2 = un

j+1/2 − ν(gRj+1 − gR

j ), j ≥ 0,

(27)



with for all j 6= 0:

gEj = gE(un

j−1/2,unj+1/2), gR

j = gR(unj−1/2,u

nj+1/2).

We concentrate on 3-point conservative schemes without loss of generality. At thediscrete level, the coupling of HEM and HRM amounts to define the quantities gE

0

and gR0 . In particular, observe from now on that the first step is conservative in the

quantities common to both systems, namely ρ, ρu and ρE, if and only if the lastthree components of gE

0 and gR0 are equal. We now describe three different strategies

for evaluating gE0 and gR

0 .

4.1.1 The flux coupling

The flux coupling strategy consists in including HEM and HRM to be coupled ina global model. The very motivation is to propose a fully conservative numericaltreatment. To that purpose, let us first recall that owing to Remark 4, HEM andHRM are thermodynamically consistent, that is pR,eq = pE. In a rather naturalway, we then consider the following HRM-like relaxation system in order to includetogether HEM and HRM in the same formalism:

∂t(ρ1z) + ∂x(ρ1zu) = λ(x)(ρ?1z

?(ρ) − ρ1z),

∂tρ + ∂x(ρu) = 0,∂t(ρu) + ∂x(ρu2 + p) = 0,∂t(ρE) + ∂x(ρE + p)u = 0,

(28)

with (x, t) ∈ R × R+,? and

λ(x) =

+∞ for x < 0,

0 for x > 0.(29)

As expected, the relaxation parameter λ is considered to be +∞ for HEM, that isρ1z = ρ?

1z?(ρ), and 0 for HRM since for the present moment we are dealing with the

convective parts only. The pressure p = pR(ρ1z, ρ, ε) is still given by (23). Definingu

n,eq−1/2 =

(ρ?1z

?(ρ), ρ, ρu, ρE)n

−1/2, we are thus led to set

gE0 = (gR

2 , gR3 , gR

4 )(un,eq−1/2,u

n1/2) and gR

0 = gR(un,eq−1/2,u

n1/2), (30)

where with classical notations, gRi i=2,3,4 denotes the last three components of gR.

4.1.2 The intermediate state coupling

We have just seen that the flux coupling is motivated by (and achieves) a conser-vation property on conservative variables ρ, ρu, ρE common to both systems HEMand HRM. Regarding the intermediate state coupling, the idea is rather to imposethe continuity of some variables of physical interest through the interface.

We first propose to impose the continuity of the common variables ρ, ρu, ρEat the interface. So, introducing the natural vectors u

n,E1/2 = (ρ, ρu, ρE)n

1/2 and



un,R−1/2 =

(ρ?1z

?(ρ), ρ, ρu, ρE)n

−1/2for converting a vector of HRM into a vector of

HEM and vice versa, it gives the following natural definition for gE0 and gR

0 :

gE0 = gE(un

−1/2,un,E1/2) and gR

0 = gR(un,R−1/2,u

n1/2). (31)

Remark 6 We remark that un,R−1/2 = u

n,eq−1/2 so that (30) and (31) only differ by the

definition of gE0 . But in the case when un

1/2 is at equilibrium, that is (ρ1z)n1/2 =

ρ?1z

?(ρn1/2), let us recall that the pressure p obeys equivalently (15) and (23) (see

Remark 4). This implies that the last three conservation laws of HRM on ρ, ρu andρE coincide with the three ones of HEM. From a numerical point of view, it is thusexpected that (gR

2 , gR3 , gR

4 )(un,eq−1/2,u

n1/2) and gE(un

−1/2,un,E1/2) are equal, or at least very

close depending on the choice of gE and gR.

Definition (31) aims at providing, whenever possible, the continuity of ρ, ρuand ρE at the interface. If such a property actually holds, the internal energy ε iscontinuous as well by (17) but not the pressure p since generally speaking un

1/2 is not

at equilibrium, that is ρ1z 6= ρ?1z

?(ρ) and therefore pE(ρ, ε) 6= pR(ρ1z, ρ, ε). Thatis the reason why we now propose to modify the initial intermediate state couplingin order to impose the continuity of the pressure p and let us say ρ and ρu. Thelatter choice is natural to achieve conservation of mass ρ and momentum ρu. Wethen define two vectors, u

n,E1/2 for HEM and u

n,R−1/2 for HRM, sharing the same ρ, ρu

and p as un1/2 and un

−1/2 respectively:

un,E1/2 = (ρ, ρu, ρE

E)n1/2 and u

n,R−1/2 =

(ρ?1z

?(ρ), ρ, ρu, ρER)n

−1/2.

This is done by inverting the pressure laws (15) and (23) with respect to ε. Usingin addition (17), straightforward calculations give:

ρEE

=1

2

(ρu)2

ρ+

p

(γ1 − 1)if ρ ≤ ρ?

1,

p

(γ1 − 1)

ρ

ρ?1

if ρ?1 < ρ < ρ?

2,

p

(γ2 − 1)if ρ ≥ ρ?

2

and

ρER

=1

2

(ρu)2

ρ+

ρp

(γ1 − 1)ρ?1z

?(ρ) + (γ2 − 1)(ρ − ρ?

1z?(ρ)

) .

Then we set:

gE0 = gE(un

−1/2,un,E1/2) and gR

0 = gR(un,R−1/2,u

n1/2). (32)

Remark 7 The state un,R−1/2 being at equilibrium, it is clear that ρE

Rn

−1/2 equals

ρEn−1/2. As an immediate consequence we have u

n,R−1/2 = u

n,R−1/2 and both flux coupling

and intermediate state coupling strategies (modified or not) only differ by gE0 (see

also Remark 6).



Remark 8 Let us also note that if un1/2 is at equilibrium, the corresponding pressure

pn1/2 may well be computed using either pE or pR so that u

n,E1/2 = u

n,E1/2 and both

intermediate state coupling strategies (modified or not) have the same gE0 . Hence

and by Remark 6, the three coupling strategies are expected to be equal in this (veryparticular) situation.

To conclude this section, let us keep in mind that in general gE0 6= gR

0 for theintermediate state coupling strategy (modified or not). Therefore, the first step ofthe whole numerical procedure is not conservative in ρ, ρu and ρE contrary to theflux coupling strategy.

4.2 The source term

This section is devoted to the numerical treatment of (26). Since SE = (0, 0, 0) andSR = (λ0

(ρ?1z

?(ρ) − ρ1z), 0, 0, 0), we naturally set

ρn+1j+1/2 = ρn+1−

j+1/2,

un+1j+1/2 = un+1−

j+1/2,

En+1j+1/2 = En+1−

j+1/2,

∀j ∈ Z, n ∈ N

since ρ, ρu and ρe do not vary in this step. It remains to solve

∂tρ1z = λ0

(ρ?1z

?(ρ) − ρ1z)

in each cell of DR. This is done exactly via the formula

ρ1n+1j+1/2

zn+1j+1/2

= ρ?1z

?(ρn+1−j+1/2

) −(ρ?1z

?(ρn+1−j+1/2

) − ρ1n+1−j+1/2

zn+1−j+1/2

)e−λ0∆t,

for all j ∈ N and n ∈ N.

4.3 The numerical schemes

We present now the numerical schemes we have tested, that is, we give severalpossible definitions for gE and gR. The first scheme is the Rusanov scheme [34],which is a very simple, but diffusive, scheme. The second scheme is a standardLagrange-Projection scheme, as proposed in [20]. In opposition to the previous one,this method is an upwinding scheme, but not a 3-point scheme (it is actually a5-point scheme and the numerical flux also depends on ν = ∆t/∆x). The thirdscheme is a nonconservative modification of the Lagrange-Projection scheme firstintroduced in [4, 8], in order to avoid spurious oscillations of pressure near contactdiscontinuities.

4.3.1 The Rusanov scheme

The Rusanov scheme [34] is a classical 3-point scheme. The associated numericalflux is

gα(u,v) =fα(u) + fα(v)

2−

λαm(u,v)

2(v − u), (33)



withλα

m(u,v) = max(maxi

(|λαi (u)|),max

i(|λα

i (v)|))

where λαi is the i-th eigenvalue of the jacobian matrix Df α, α = E,R. Under the

classical CFL condition

ν maxj

(maxi

|λαi (un

j+1/2)|) ≤ C <1

2,

the Rusanov scheme is positive for the density and stable (C is called the Courantnumber). This scheme is very simple but it is very diffusive, which is due to thecentral, but stable, discretization (opposed to an upwind discretization). When con-sidering for instance stationary contact discontinuities (null velocity, uniform pres-sure and non uniform density), the Rusanov scheme cannot preserve such profiles,the density is diffused, contrarily to usual upwind schemes like Godunov method orLagrange-Projection schemes. We will see the consequences of this drawback in thenumerical results.

4.3.2 The Lagrange-Projection scheme

Lagrange-Projection schemes are based on an operator splitting consisting in solv-ing first the Euler equations in pseudo-Lagrangian coordinates (tn → tn+1/2) andthen the advection part of the equations (tn+1/2 → tn+1−), which may appear as aprojection procedure on the fixed mesh. The exact Godunov scheme associated withthis splitting is described in [24]. We prefer here to use an approximate Godunovresolution of the Lagrange part. This Lagrangian scheme can be interpreted as anacoustic scheme (see [21]) or a relaxation scheme (see [19]). The Lagrange step ofthe scheme for HRM is

ρ1n+1/2j+1/2 z

n+1/2j+1/2 = ρ

n+1/2j+1/2

ρ1nj+1/2z

nj+1/2

ρnj+1/2

,

ρn+1/2j+1/2 =

ρnj+1/2

1 + ν(un

j+1 − unj

) ,

un+1/2j+1/2 = un

j+1/2 −ν

ρnj+1/2

(pn

j+1 − pnj

),

En+1/2j+1/2 = En

j+1/2 −ν

ρnj+1/2

(pn

j+1unj+1 − pn

j unj

)

with un

j = 12(ρc)n

j(pn

j−1/2 − pnj+1/2) + 1

2(unj−1/2 + un

j+1/2),

pnj = 1

2(pnj−1/2 + pn

j+1/2) +(ρc)n

j

2 (unj−1/2 − un

j+1/2),

the approximate local Lagrangian sound speed (ρc)nj above being given by

(ρc)nj =

√max(ρn

j−1/2cnj−1/2

2, ρnj+1/2c

nj+1/2

2)min(ρnj−1/2, ρ

nj+1/2)

where cnj+1/2 is the sound speed in cell j + 1/2 at time step n: cn

j+1/2 = cR(unj+1/2).

Quantities pnj+1/2 are to be understood in the same sense,

pnj+1/2 = pR

(ρ1

nj+1/2z

nj+1/2, ρ

nj+1/2, E

nj+1/2 − 1/2

(un

j+1/2

)2)



for all j under consideration.Then, the projection step enables to compute

(ρ1

n+1−j+1/2z

n+1−j+1/2, ρ

n+1−j+1/2, ρ

n+1−j+1/2u

n+1−j+1/2, ρ

n+1−j+1/2E

n+1−j+1/2

).

The whole Lagrange-Projection time step for HRM results in the following flux:

gRj =

ρ1nj zn

j unj

ρnj un

j

ρnj un

j unj + pn

j

ρnj En

j unj + pn

j unj

(see above for the fluxes unj and pn

j ) with

if unj ≥ 0,

ρ1nj zn

j = ρ1n+1/2j−1/2 z

n+1/2j−1/2 ,

ρnj = ρ

n+1/2j−1/2 ,

unj = u

n+1/2j−1/2

,

Enj = E

n+1/2j−1/2 ,

and if unj < 0,

ρ1nj zn

j = ρ1n+1/2j+1/2 z

n+1/2j+1 ,

ρnj = ρ

n+1/2j+1/2 ,

unj = u

n+1/2j+1/2

,

Enj = E

n+1/2j+1/2 .

The flux gEj for HEM follows the same definition and is composed of the last three

components of gRj .

Note that because of the definition of ρ1nj zn

j , ρnj , un

j and Enj , this is a 4-point

flux. This leads to define four vectors instead of two for the numerical coupling:u

n,R−3/2 and u

n,R−1/2 for HRM, and u

n,E1/2 and u

n,E3/2 for HEM.

This scheme has the drawback of creating spurious oscillations near contact dis-continuities when dealing with a non-ideal gas such as in HEM and HRM. Forthis reason we propose a slight modification of the Lagrange-Projection scheme toavoid this phenomenon. The modification is based on the fact that the oscilla-tions come up in the projection procedure, the Lagrange step being oscillation-free whatever the pressure law. Thus we propose to project the quantities ρ1z,ρ, ρu, and p instead of ρE,which implies a maximum principle on p in the pro-jection step. This new scheme is oscillation-free near contact discontinuities buthas the drawback of being nonconservative in total energy. Nevertheless in thefollowing numerical test-cases where shocks are not too strong, we do not ob-serve important artefacts. This scheme will be called “nonconservative Lagrange-Projection scheme” and noted L-PP in the following. Let us briefly give the schemefor HRM. The only difference concerns the quantity ρE. It is updated in twosteps. In the Lagrange step, the same formulas as above are used, and we com-

pute pn+1/2j+1/2 =

(ρ1

n+1/2j+1/2 z

n+1/2j+1/2 , ρ

n+1/2j+1/2 , E

n+1/2j+1/2 − 1/2

(u

n+1/2j+1/2

)2)

. Then, the un-

knowns (ρ1z, ρ, ρu) are updated in the same way as with the conservative Lagrange-Projection scheme, but not E. We here compute first pn+1−

j+1/2as

pn+1−j+1/2

= pn+1/2j+1/2

− ν(un

j+1(pnj+1 − p

n+1/2j+1/2

) − unj (pn

j − pn+1/2j+1/2

))

withif un

j ≥ 0, pnj = p

n+1/2j−1/2 , and if un

j < 0, pnj = p

n+1/2j+1/2 .



Finally, we compute En+1−j+1/2 by inverting the pressure law, finding En+1−

j+1/2 such that

pR

(ρ1

n+1−j+1/2z

n+1−j+1/2, ρ

n+1−j+1/2, E

n+1−j+1/2 − 1/2

(un+1−

j+1/2

)2)

= pn+1−j+1/2.

The nonconservative scheme for HEM is a straightforward transposition of thisone.

5 Numerical experiments

In this section, several numerical tests are presented in order to illustrate the differentbehaviors obtained at the coupling interface, according to the numerical scheme andthe numerical coupling we use. In all the experiments, the computational spacedomain is [−1/2; 1/2], the associated mesh is composed of 500 cells and the Courantnumber C equals 0.4. We will always consider Riemann initial data given by

uE(0, x) = ul, −1/2 ≤ x < 0,uR(0, x) = ur, 0 < x ≤ 1/2,

(34)

where ul ∈ ΩE and ur ∈ ΩR are two constant states. Actually, the initial datawill be given with respect to the variable v = (ρ, u, p). More precisely, ul =(ρl, ρlul, ρ(εE(ρl, pl) + (ul)

2/2)) where εE(ρ, p) is given by inverting (15) and ur =((ρ1)rzr, ρr, ρrur, ρ(εR((ρ1)rzr, ρr, pr) + (ur)

2/2)) where

(ρ1)r =ρr

zr + γ1−1γ2−1(1 − zr)

, (35)

and εR(ρ1z, ρ, p) is provided by (23). Besides, the initial data for the HRM part willbe given with respect to the variable (c,v), where c is the mass fraction of vapor andis plotted in each figure. It is defined by c = ρ?

1z?(ρ)/ρ for x < 0 and by c = ρ1z/ρ

for x > 0, and it naturally lies between 0 and 1.The specific heat Cv is equal to 1, and the values of the adiabatic coefficients are

γ1 = 1.6 and γ2 = 1.4, which leads to ρ?1 ≈ 0.613132 and ρ?

2 ≈ 0.919699, using (13)(below, these values will be plotted on each graph of the density variable).

The first five tests illustrate the ability of the different coupling methods toprovide the continuity at the coupling interface of (ρu, ρu2 + p, u(ρE + p)) for theflux coupling, of (ρ, ρu, ρE) for the intermediate state coupling with the conservativevariable and of (ρ, u, p) for the modified intermediate state coupling. In order tomake the results we present as clear as possible, the source term of HRM is not takeninto account in these five experiments. In the sixth one, the numerical convergenceof the coupling model between HEM and HRM toward a global HEM, letting λ0

increase, is investigated.



5.1 A simple test in phase 2

We begin these numerical tests with a very simple case. The data of this test are

c ρ u p

vl − 2 0 1

(c,v)r 0 1.5 0 2

tm = 0.2, λ0 = 0,

(36)

where tm represents the time at which the solutions are plotted. One may remarkthat the HRM part of the domain is initially at equilibrium (zr = z?(ρr)), so thatHRM becomes equivalent to HEM. Moreover, let us emphasize that this solution istotally involved in phase 2 (that is ρ > ρ?

2). Therefore, the global solution of thecoupling problem corresponds to the one provided by a unique system HEM, definedeverywhere in the domain of computation. Due to these particular properties of theinitial data, the three numerical coupling methods must give the same results (seealso Rem. 6 and 8), which may be seen in Figs. 1, 2 and 3. One may also checkthat the Rusanov scheme is more diffusive than the Lagrange-Projection schemes.

1

1.2

1.4

1.6

1.8

2

-0.4 -0.2 0 0.2 0.4

Pressure

L-PL-PPRus

-0.3

-0.25

-0.2

-0.15

-0.1

-0.05

0

-0.4 -0.2 0 0.2 0.4

Velocity

L-PL-PPRus

0.5

1

1.5

2

2.5

3

-0.4 -0.2 0 0.2 0.4

Density

L-PL-PPRus

-1

-0.5

0

0.5

1

-0.4 -0.2 0 0.2 0.4

Fraction of vapor

L-PL-PPRus

Figure 1: A test in phase 2 (36): Intermediate state coupling (conservative vari-able). Lagrange-Projection (L-P), nonconservative Lagrange-Projection (L-PP) andRusanov scheme (Rus).



1

1.2

1.4

1.6

1.8

2

-0.4 -0.2 0 0.2 0.4

Pressure

L-PL-PPRus

-0.3

-0.25

-0.2

-0.15

-0.1

-0.05

0

-0.4 -0.2 0 0.2 0.4

Velocity

L-PL-PPRus

0.5

1

1.5

2

2.5

3

-0.4 -0.2 0 0.2 0.4

Density

L-PL-PPRus

-1

-0.5

0

0.5

1

-0.4 -0.2 0 0.2 0.4

Fraction of vapor

L-PL-PPRus

Figure 2: A test in phase 2 (36): Intermediate state coupling (nonconservativevariable). Lagrange-Projection (L-P), nonconservative Lagrange-Projection (L-PP)and Rusanov scheme (Rus).



1

1.2

1.4

1.6

1.8

2

-0.4 -0.2 0 0.2 0.4

Pressure

L-PL-PPRus

-0.3

-0.25

-0.2

-0.15

-0.1

-0.05

0

-0.4 -0.2 0 0.2 0.4

Velocity

L-PL-PPRus

0.5

1

1.5

2

2.5

3

-0.4 -0.2 0 0.2 0.4

Density

L-PL-PPRus

-1

-0.5

0

0.5

1

-0.4 -0.2 0 0.2 0.4

Fraction of vapor

L-PL-PPRus

Figure 3: A test in phase 2 (36): Flux coupling. Lagrange-Projection (L-P), non-conservative Lagrange-Projection (L-PP) and Rusanov scheme (Rus).



5.2 Uniform velocity and pressure at equilibrium

In this case, the right part of the domain is still at equilibrium, but the left partand the right part belong to different phases since ρl > ρ?

2 and ρr < ρ?1. Since

zr = z?(ρr), this test enters in the frame described in Rem. 6 and 8. The numericalresults of the different coupling methods are thus expected to be very close. Thedata associated with this test are

c ρ u p

vl − 2 1 1

(c,v)r 1 0.5 1 1

tm = 0.15, λ0 = 0,

(37)

so that the solution is composed by a contact discontinuity, moving to the right,that is in the HRM part of the domain. Since the contact discontinuity separatesthe two phases, the adiabatic coefficients are different on both sides of the waveand it is well-known that, in such case, any standard conservative scheme providesspurious oscillations on the profiles of u and p (see for instance [1]). Basically, only

0.97

0.975

0.98

0.985

0.99

0.995

1

1.005

-0.4 -0.2 0 0.2 0.4

Pressure

L-PL-PPRus

0.97 0.975

0.98

0.985 0.99

0.995 1

1.005 1.01

1.015

-0.4 -0.2 0 0.2 0.4

Velocity

L-PL-PPRus

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

-0.4 -0.2 0 0.2 0.4

Density

L-PL-PPRus

0

0.2

0.4

0.6

0.8

1

-0.4 -0.2 0 0.2 0.4

Fraction of vapor

L-PL-PPRus

Figure 4: Uniform velocity and pressure at equilibrium (37): Intermediate state cou-pling (conservative variable). Lagrange-Projection (L-P), nonconservative Lagrange-Projection (L-PP) and Rusanov scheme (Rus).

the nonconservative Lagrange-Projection scheme can maintain the velocity and thepressure constant, that is exactly the reason why this scheme has been introduced(see also [4, 8]). Nevertheless, as mentioned above, the intermediate state coupling



0.97

0.975

0.98

0.985

0.99

0.995

1

-0.4 -0.2 0 0.2 0.4

Pressure

L-PL-PPRus

0.97 0.975

0.98

0.985 0.99

0.995 1

1.005 1.01

1.015

-0.4 -0.2 0 0.2 0.4

Velocity

L-PL-PPRus

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

-0.4 -0.2 0 0.2 0.4

Density

L-PL-PPRus

0

0.2

0.4

0.6

0.8

1

-0.4 -0.2 0 0.2 0.4

Fraction of vapor

L-PL-PPRus

Figure 5: Uniform velocity and pressure at equilibrium (37): Intermediate statecoupling (nonconservative variable). Lagrange-Projection (L-P), nonconservativeLagrange-Projection (L-PP) and Rusanov scheme (Rus).



0.97

0.975

0.98

0.985

0.99

0.995

1

-0.4 -0.2 0 0.2 0.4

Pressure

L-PL-PPRus

0.965 0.97

0.975 0.98

0.985 0.99

0.995 1

1.005 1.01

1.015

-0.4 -0.2 0 0.2 0.4

Velocity

L-PL-PPRus

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

-0.4 -0.2 0 0.2 0.4

Density

L-PL-PPRus

0

0.2

0.4

0.6

0.8

1

-0.4 -0.2 0 0.2 0.4

Fraction of vapor

L-PL-PPRus

Figure 6: Uniform velocity and pressure at equilibrium (37): Flux coupling.Lagrange-Projection (L-P), nonconservative Lagrange-Projection (L-PP) and Ru-sanov scheme (Rus).



with the conservative variable modifies the pressure through the coupling interface.Consequently, variations of pressure and velocity can be observed in Fig. 4 even forthe nonconservative Lagrange-Projection scheme L-PP.

5.3 Uniform conservative variables out of equilibrium

We now focus on initial data out of equilibrium, that is zr 6= z?(ρr). The initial dataof this test case is

c ρ u p

vl − 2 1 1

(c,v)r 1 2 1 1.5

tm = 0.15, λ0 = 0.

(38)

One may check by a simple calculation that (ρl, (ρu)l, (ρE)l) = (ρr, (ρu)r, (ρE)r)and z?(ρl) 6= zr. In both models, the evolution in time of the conservative variables

2.6

2.8

3

3.2

3.4

3.6

3.8

4

-0.4 -0.2 0 0.2 0.4

Energy

L-PL-PPRus

1.5 1.55

1.6 1.65

1.7 1.75

1.8 1.85

1.9 1.95

2 2.05

-0.4 -0.2 0 0.2 0.4

Impulsion

L-PL-PPRus

0.6 0.8

1

1.2 1.4 1.6 1.8

2 2.2 2.4

-0.4 -0.2 0 0.2 0.4

Density

L-PL-PPRus

0

0.2

0.4

0.6

0.8

1

-0.4 -0.2 0 0.2 0.4

Fraction of vapor

L-PL-PPRus

Figure 7: Uniform conservative variables out of equilibrium (38): Intermediatestate coupling (conservative variable). Lagrange-Projection (L-P), nonconservativeLagrange-Projection (L-PP) and Rusanov scheme (Rus).

ρ, ρu and ρE is partially due to the spatial variation of the pressure ∂xp. Here,the pressure is initially discontinuous at the coupling interface (because z?(ρl) 6= zr)and thus ρ, ρu and ρE cannot remain constant for t > 0, even when using theintermediate state coupling with the conservative variable, see Fig. 7. Besides, sincethe velocity is positive, the discontinuity moves to the right and after some time



1

1.1

1.2

1.3

1.4

1.5

-0.4 -0.2 0 0.2 0.4

Pressure

L-PL-PPRus

0.86

0.88

0.9

0.92

0.94

0.96

0.98

1

-0.4 -0.2 0 0.2 0.4

Velocity

L-PL-PPRus

0.6 0.8

1

1.2 1.4 1.6 1.8

2 2.2 2.4

-0.4 -0.2 0 0.2 0.4

Density

L-PL-PPRus

0

0.2

0.4

0.6

0.8

1

-0.4 -0.2 0 0.2 0.4

Fraction of vapor

L-PL-PPRus

Figure 8: Uniform conservative variables out of equilibrium (38): Intermediate statecoupling (nonconservative variable). Lagrange-Projection (L-P), nonconservativeLagrange-Projection (L-PP) and Rusanov scheme (Rus).



3.4

3.6

3.8

4

4.2

4.4

4.6

4.8

5

-0.4 -0.2 0 0.2 0.4

Energy flux

L-PL-PPRus

2.5 2.6 2.7 2.8 2.9

3 3.1 3.2 3.3 3.4 3.5

-0.4 -0.2 0 0.2 0.4

Impulsion flux

L-PL-PPRus

1.5 1.55

1.6 1.65

1.7 1.75

1.8 1.85

1.9 1.95

2 2.05

-0.4 -0.2 0 0.2 0.4

Density flux

L-PL-PPRus

0

0.2

0.4

0.6

0.8

1

-0.4 -0.2 0 0.2 0.4

Fraction of vapor

L-PL-PPRus

Figure 9: Uniform conservative variables out of equilibrium (38): Flux coupling.Lagrange-Projection (L-P), nonconservative Lagrange-Projection (L-PP) and Ru-sanov scheme (Rus).



steps, the cell [x0, x1] is at equilibrium (i.e. zn1/2 = z?(ρn

1/2)), leading to a continuouspressure law through the coupling interface x = 0. As a result, all the variables arecontinuous through this interface (see Figs. 7, 8 and 9, where we have plotted thevariable corresponding to the coupling method we used).

5.4 Uniform primitive variables out of equilibrium

In this case, the HRM part is still out of equilibrium:

c ρ u p

vl − 1 −0.5 1

(c,v)r 1 1 −0.5 1

tm = 0.2, λ0 = 0.

(39)

Let us note that the density, the velocity and the pressure (that are the primitivevariables) are the same on both sides of the coupling interface. In fact, if the state(c,v)r was at equilibrium, the primitive variables would be preserved constant byall the numerical methods, as we saw in the results of Section 5.2. But, for thistest, the equation of state changes at the coupling interface (since zr 6= z?(ρr)) andsolutions with complex structures can be developed.

0.75

0.8

0.85

0.9

0.95

1

-0.4 -0.2 0 0.2 0.4

Pressure

L-PL-PPRus

-0.55

-0.5

-0.45

-0.4

-0.35

-0.3

-0.4 -0.2 0 0.2 0.4

Velocity

L-PL-PPRus

0.6

0.7

0.8

0.9

1

1.1

1.2

1.3

-0.4 -0.2 0 0.2 0.4

Density

L-PL-PPRus

0

0.2

0.4

0.6

0.8

1

-0.4 -0.2 0 0.2 0.4

Fraction of vapor

L-PL-PPRus

Figure 10: Uniform velocity and pressure out of equilibrium (39): Intermediatestate coupling (conservative variable). Lagrange-Projection (L-P), nonconservativeLagrange-Projection (L-PP) and Rusanov scheme (Rus).



0.995 1

1.005 1.01

1.015 1.02

1.025 1.03

1.035 1.04

1.045

-0.4 -0.2 0 0.2 0.4

Pressure

L-PL-PPRus

-0.515-0.51

-0.505

-0.5-0.495-0.49

-0.485

-0.48-0.475-0.47

-0.4 -0.2 0 0.2 0.4

Velocity

L-PL-PPRus

0.6 0.65

0.7

0.75 0.8

0.85 0.9

0.95 1

1.05

-0.4 -0.2 0 0.2 0.4

Density

L-PL-PPRus

0

0.2

0.4

0.6

0.8

1

-0.4 -0.2 0 0.2 0.4

Fraction of vapor

L-PL-PPRus

Figure 11: Uniform velocity and pressure out of equilibrium (39): Intermediate statecoupling (nonconservative variable). Lagrange-Projection (L-P), nonconservativeLagrange-Projection (L-PP) and Rusanov scheme (Rus).



0.7

0.75

0.8

0.85

0.9

0.95

1

-0.4 -0.2 0 0.2 0.4

Pressure

L-PL-PPRus

-0.7

-0.65

-0.6

-0.55

-0.5

-0.45

-0.4

-0.4 -0.2 0 0.2 0.4

Velocity

L-PL-PPRus

0.6

0.7

0.8

0.9

1

1.1

1.2

1.3

1.4

-0.4 -0.2 0 0.2 0.4

Density

L-PL-PPRus

0

0.2

0.4

0.6

0.8

1

-0.4 -0.2 0 0.2 0.4

Fraction of vapor

L-PL-PPRus

Figure 12: Uniform velocity and pressure out of equilibrium (39): Flux coupling.Lagrange-Projection (L-P), nonconservative Lagrange-Projection (L-PP) and Ru-sanov scheme (Rus).



Only the two Lagrange-Projection schemes with the intermediate state couplingbased on the nonconservative variable (ρ, ρu, p) preserve the primitive variables con-stant (see Fig. 11), all other methods introduce acoustic waves in the solution (Figs.10, 11 and 12), notice the mixture zone in Fig. 10). More precisely, since the equa-tion of state is different on both sides of the coupling interface, the flux couplingand the intermediate state coupling based on the conservative variable provide dif-ferent values of the pressure from the one side to the other side and thus, introduceacoustic waves.

The reason why the Rusanov scheme with the modified intermediate state cou-pling introduces acoustic waves is different. As we mentioned, the Rusanov schemeis very diffusive. Therefore, the coupling interface is slightly “diffused” by the Ru-sanov scheme, in the sense that, in the first cell at the right of the coupling interface,c no longer exactly equals either 0 or 1 but lies strictly in (0, 1) (see the shape of thefraction of vapor in Fig. 11). The direct consequence is a modification of the pres-sure in the cell at the right of the coupling interface and thus, the pressure cannotbe left constant.

5.5 Shock tube with occurrence of phase 1

In this test, the initial condition is at equilibrium:

c ρ u p

vl − 1 −2 1

(c,v)r 0 1 1 1

tm = 0.1, λ0 = 0,

(40)

but, whereas the initial condition is in phase 2, an intermediate zone with ρ < ρ?1

appears for t > 0, overlapping the coupling interface. Therefore, c = 1 in the HEMpart of the intermediate zone but since the velocity at x = 0 is negative, the fractionc stays equal to 0 in the HRM part of the intermediate zone. As in the previous case,the discontinuity of z leads to a discontinuity of the pressure law. Therefore, onlythe coupling method based on the nonconservative variable provides a continuouspressure through x = 0, see Fig. 14. In Fig. 15 are plotted the results obtainedby the flux coupling. The variables ρu, ρu2 + p and u(ρE + p) are represented andone may see that, as expected, they are constant through the coupling interface.However, in Fig. 13, the variables ρ, ρu and ρE are discontinuous at x = 0, thoughthese results correspond to the intermediate state coupling with the conservativevariable. Indeed, we have seen in Sec. 5.3 (see also Fig. 7) that if the HRM part isout of equilibrium, the continuity of the conservative variables cannot be achieved(in test of Sec. 5.3, this leads to the appearance of acoustic waves). In the presentcase, since c, and thus z, is discontinuous at x = 0 (for all t > 0), the conservativevariables are discontinuous at the coupling interface.

5.6 Convergence of HRM towards equilibrium

This test is based on a uniform velocity and a uniform pressure. Moreover, the HRMpart is out of equilibrium (this test is similar to the test of Sec. 5.4, except for the



0 0.5

1

1.5 2

2.5 3

3.5 4

4.5

-0.4 -0.2 0 0.2 0.4

Energy

L-PL-PPRus

-2

-1.5

-1

-0.5

0

0.5

1

-0.4 -0.2 0 0.2 0.4

Impulsion

L-PL-PPRus

0.1 0.2 0.3

0.4 0.5 0.6 0.7

0.8 0.9

1

-0.4 -0.2 0 0.2 0.4

Density

L-PL-PPRus

0

0.2

0.4

0.6

0.8

1

-0.4 -0.2 0 0.2 0.4

Fraction of vapor

L-PL-PPRus

Figure 13: Shock tube with occurrence of phase 1 (40): Intermediate state cou-pling (conservative variable). Lagrange-Projection (L-P), nonconservative Lagrange-Projection (L-PP) and Rusanov scheme (Rus).



0.1 0.2 0.3

0.4 0.5 0.6 0.7

0.8 0.9

1

-0.4 -0.2 0 0.2 0.4

Pressure

L-PL-PPRus

-2

-1.5

-1

-0.5

0

0.5

1

-0.4 -0.2 0 0.2 0.4

Velocity

L-PL-PPRus

0.1 0.2 0.3

0.4 0.5 0.6 0.7

0.8 0.9

1

-0.4 -0.2 0 0.2 0.4

Density

L-PL-PPRus

0

0.2

0.4

0.6

0.8

1

-0.4 -0.2 0 0.2 0.4

Fraction of vapor

L-PL-PPRus

Figure 14: Shock tube with occurrence of phase 1 (40): Intermediate state coupling(nonconservative variable). Lagrange-Projection (L-P), nonconservative Lagrange-Projection (L-PP) and Rusanov scheme (Rus).



-12

-10

-8

-6

-4

-2

0

2

4

-0.4 -0.2 0 0.2 0.4

Energy flux

L-PL-PPRus

0 0.5

1 1.5

2 2.5

3 3.5

4 4.5

5

-0.4 -0.2 0 0.2 0.4

Impulsion flux

L-PL-PPRus

-2

-1.5

-1

-0.5

0

0.5

1

-0.4 -0.2 0 0.2 0.4

Density flux

L-PL-PPRus

0

0.2

0.4

0.6

0.8

1

-0.4 -0.2 0 0.2 0.4

Fraction of vapor

L-PL-PPRus

Figure 15: Shock tube with occurrence of phase 1 (40): Flux coupling. Lagrange-Projection (L-P), nonconservative Lagrange-Projection (L-PP) and Rusanov scheme(Rus).



density). We use different values of the relaxation parameter λ0 in HRM:

c ρ u p

vl − 1 −0.5 1

(c,v)r 1 2 −0.5 1

tm = 0.2, λ0 = 0, 10, 100.

(41)

As we mentioned in Rem. 4, HRM formally tends to HEM when λ0 → +∞. This isthe result we want to test, adding the difficulty of the numerical coupling at x = 0.But, before commenting the numerical results, let us study the solution expected forlarge λ0. One may think naively that the limit solution would be the HEM solutionassociated with the initial data

ρ u p

vl 1 −0.5 1

vr 2 −0.5 1

(42)

But, for the coupling problem, zr tends to z?(ρr). As a consequence, the pressure ismodified in the HRM part since pR depends on ρ1z and, as noted before, (ρ1z)r 6=ρ?1z

?(ρr). On the contrary, the pressure of the solution of HEM with data (42)remains constant (in time and space). Then, what is the limit solution? In fact, itis the HEM solution with the initial data

ρ u p

vl 1 −0.5 1

vr 2 −0.5 2/3

(43)

since, using Rem. 4 and (35), we have

pr = pR,eq(ρr, εR((ρ1)rz

?(ρr), ρr, pr)),

= pR,eq(2, εR(0, 2, 1)),

= 2/3.

In Figs. 16, 17 and 18 are plotted the results for several λ0, using the Lagrange-Projection scheme. One may see that, as expected, all the coupling methods tendsto the solution based on HEM with data (43), whatever the numerical scheme andthe coupling method.

6 Conclusion

We have shown that the coupling problem (1) can be numerically solved, once it iscompleted with an interface model that restores the continuity of physics. It must benoted that there is a multiple choice of interface models that can apply, dependingon the physics that is under study. Moreover, in each case, we are able to givea numerical treatment that will verify the constraints imposed by the interfacialmodel, at least in a weak sense. The application of the interface coupling developed



0.65

0.7

0.75

0.8

0.85

0.9

0.95

1

-0.4 -0.2 0 0.2 0.4

Pressure

HEM0

10100

-0.6

-0.55

-0.5

-0.45

-0.4

-0.35

-0.4 -0.2 0 0.2 0.4

Velocity

HEM0

10100

0.6 0.8

1 1.2 1.4 1.6 1.8

2 2.2 2.4 2.6

-0.4 -0.2 0 0.2 0.4

Density

HEM0

10100

0

0.2

0.4

0.6

0.8

1

-0.4 -0.2 0 0.2 0.4

Fraction of vapor

HEM0

10100

Figure 16: Convergence towards equilibrium (41): Intermediate state coupling (con-servative variable). Lagrange-Projection scheme for different λ0.



0.65

0.7

0.75

0.8

0.85

0.9

0.95

1

-0.4 -0.2 0 0.2 0.4

Pressure

HEM0

10100

-0.5

-0.48

-0.46

-0.44

-0.42

-0.4

-0.38

-0.36

-0.4 -0.2 0 0.2 0.4

Velocity

HEM0

10100

0.6 0.8

1 1.2 1.4 1.6 1.8

2 2.2 2.4 2.6

-0.4 -0.2 0 0.2 0.4

Density

HEM0

10100

0

0.2

0.4

0.6

0.8

1

-0.4 -0.2 0 0.2 0.4

Fraction of vapor

HEM0

10100

Figure 17: Convergence towards equilibrium (41): Intermediate state coupling (non-conservative variable). Lagrange-Projection scheme for different λ0.



0.4

0.5

0.6

0.7

0.8

0.9

1

-0.4 -0.2 0 0.2 0.4

Pressure

HEM0

10100

-0.8-0.75

-0.7

-0.65-0.6

-0.55-0.5

-0.45-0.4

-0.35

-0.4 -0.2 0 0.2 0.4

Velocity

HEM0

10100

0.6 0.8

1 1.2 1.4 1.6 1.8

2 2.2 2.4 2.6

-0.4 -0.2 0 0.2 0.4

Density

HEM0

10100

0

0.2

0.4

0.6

0.8

1

-0.4 -0.2 0 0.2 0.4

Fraction of vapor

HEM0

10100

Figure 18: Convergence towards equilibrium (41): Flux coupling. Lagrange-Projection scheme for different λ0.



here to other types of numerical methods is not the aim of this paper but seems apriori possible.

At last, let us mention that the numerical coupling techniques presented here arebeing developed in other contexts: gas dynamics (cf. [4], [5]) and Lagrangian modelsof gas dynamics (in [7]) and that the theory of the interface coupling is under study,based on the pioneering works [25] and [26]. In particular, the analysis of solutionsto Riemann problems for the coupling of systems of gas dynamics is developed in[14] and the coupling of HEM and HRM models is analyzed from a theoretical pointof view in [12].

Acknowledgements

This work was partially supported by the NEPTUNE project, funded by CEA, EDF,IRSN and AREVA-NP.

References

[1] R. Abgrall. How to prevent pressure oscillations in multicomponent flow cal-culations: a quasi-conservative approach. J. Comput. Phys., 125(1):150–160,1996.

[2] R. Abgrall, S. Karni. Computation of compressible multifluids. J. Comput.Phys., 169:594–623, 2001.

[3] G. Allaire, G. Faccanoni, S. Kokh. A Strictly Hyperbolic Equilibrium Phase Tran-sition Model, C. R. Acad. Sci. Paris Serie I, 344, pp.135-140, 2007.

[4] A. Ambroso, Ch. Chalons, F. Coquel, E. Godlewski, J.-M. Herard, F. Lagoutiere,P.-A. Raviart, N. Seguin. The coupling of multiphase flow models. Proceedingsof Nureth-11, Avignon, France, 2005.

[5] A. Ambroso, Ch. Chalons, F. Coquel, E. Godlewski, F. Lagoutiere, P.-A. Raviart,N. Seguin. Couplage de deux systemes de la dynamique des gaz. Proceedings ofCongres Francais de Mecanique, Troyes, France, 2005.

[6] A. Ambroso, Ch. Chalons, F. Coquel, E. Godlewski, F. Lagoutiere, P.-A. Raviart,N. Seguin. Homogeneous models with phase transition: coupling by Finite Vol-ume methods. Proceedings of Finite volumes for complex applications IV (Mar-rakech, 2005), Hermes Science plubisher, 483–392, 2005.

[7] A. Ambroso, Ch. Chalons, F. Coquel, E. Godlewski, F. Lagoutiere, P.-A. Raviart,N. Seguin. Extension of interface coupling to general Lagrangian systems. Pro-ceedings of ENUMATH, Santiago de Compostela, Spain, 2005.

[8] T. Barberon, P. Helluy, S. Rouy. Practical computation of axisymmetri-cal multifluid flows. Int. J. Finite Volumes, vol. 1, num. 1, 1–34, 2004(http://www.latp.univ-mrs.fr/IJFV).



[9] F. Bouchut. Nonlinear stability of finite volume methods for hyperbolic con-servation laws, and well-balanced schemes for source. Frontiers in MathematicsSeries, Birkhauser, 2004.

[10] Z. Bilicki, J. Kestin. Physical aspects of the relaxation model in two-phase flow.Proc. R. Soc. Lond., A428(1875):379–397, 1990.

[11] B. Boutin. Couplage de systemes de lois de conservation scalaires par uneregularisation de Dafermos. Master dissertation, Paris, 2005.

[12] F. Caetano. Sur certains problemes de linearisation et de couplage pour lessystemes hyperboliques non lineaires. PhD Thesis, Universite Pierre et MarieCurie-Paris6, 2006.

[13] F. Caro. Modelisation et simulation numerique des transitions de phase liquidevapeur. PhD thesis, CMAP, Ecole Polytechnique, 2004.

[14] Ch. Chalons, P.-A. Raviart and N. Seguin. The interface coupling of the gasdynamics equations. In preparation.

[15] F. Caro, F. Coquel, D. Jamet, and S. Kokh. DINMOD: A diffuse interfacemodel for two-phase flows modelling. In S. Cordier, T. Goudon, M. Gutnic,and E. Sonnendrucker, editors, Numerical methods for hyperbolic and kineticproblems, IRMA lecture in mathematics and theoretical physics (Proceedings ofthe CEMRACS 2003), 2005.

[16] F. Caro, F. Coquel, D. Jamet, and S. Kokh. Phase change simulation for isother-mal compressible two-phase flows. In AIAA Comp. Fluid Dynamics, numberAIAA-2005-4697, 2005.

[17] F. Caro, F. Coquel, D. Jamet, and S. Kokh. A Simple Finite-Volume Methodfor Compressible Isothermal Two-Phase Flows Simulation. Int. J. of Finite Vol.,vol. 3, num. 1, 2006 (http://www.latp.univ-mrs.fr/IJFV).

[18] Ch. Chalons, F. Coquel. Navier-Stokes equations with several independent pres-sure laws and explicit predictor-corrector schemes. Numer. Math., 101(3):451-478, 2005.

[19] F. Coquel, E. Godlewski, B. Perthame, A. In, M. Rascle. Some new Godunovand relaxation methods for two-phase flow problems. Godunov methods (Oxford,1999), 179–188, Kluwer/Plenum, New York, 2001.

[20] B. Despres. Symetrisation en variable de Lagrange pour la mecanique desmilieux continus et schemas numeriques. Matapli, 72:45–61, 2003.

[21] B. Despres. Invariance properties of Lagrangian systems of conservationlaws, approximate Riemann solvers and the entropy condition. Numer. Math.,89(1):99–134, 2001.



[22] P. Downar-Zapolski, Z. Bilicki, L. Bolle and J. Franco. The non equili-brium model for one-dimensionnal flashing liquid flow. Int. J. Multiphase Flow,22(3):473–483, 1996.

[23] G. Faccanoni. Modelisation fine d’ecoulements diphasiques : contribution al’etude de la crise d’ebullition. Ph.D. Thesis, Ecole Polytechnique, in prepa-ration.

[24] E. Godlewski and P.-A. Raviart. Numerical approximation of hyperbolic sys-tems of conservation laws. Applied Mathematical Sciences, 118. Springer-Verlag,New York, 1996.

[25] E. Godlewski and P.-A. Raviart. The numerical interface coupling of nonlin-ear hyperbolic systems of conservation laws. I. The scalar case. Numer. Math.,97(1):81–130, 2004.

[26] E. Godlewski, K.-C. Le Thanh and P.-A. Raviart. The numerical interfacecoupling of nonlinear hyperbolic systems of conservation laws: II. The case ofsystems. RAIRO Model. Math. Anal. Numer, 39(4):649–692, 2005.

[27] J.-M. Herard. Schemes to couple flows between free and porous medium. AIAApaper 2005-4861, 2005.

[28] J.-M. Herard, O. Hurisse. Coupling two and one-dimensional unsteady Eulerequations through a thin interface. Computer and Fluids, vol. 36, num. 4, 651-666, 2007.

[29] J.-M. Herard, O. Hurisse. A method to couple two-phase flow models. EDFInternal Report HI-81/06/001/A, France.

[30] S. Jaouen. Etude mathematique et numerique de stabilite pour des modeleshydrodynamiques avec transition de phase. PhD Thesis, Universite Pierre etMarie Curie-Paris6, 2001.

[31] Ph. G. LeFloch. Hyperbolic Systems of Conservation Laws: The theory of clas-sical and nonclassical shock waves. E.T.H. Lecture Notes Series, Birkhauser,2002.

[32] R. P. Fedkiw, T. Aslam, B. Merriman, S. Osher. A non-oscillatory Eulerian ap-proach to interfaces in multimaterial flows (the ghost fluid method). J. Comput.Phys., 152(2):457–492, 1999.

[33] J. Piraux, B. Lombard. A new interface method for hyperbolic problems withdiscontinuous coefficients: one-dimensional acoustic example. J. Comp. Phys.,168(1):227–248, 2001.

[34] V. V. Rusanov. The calculation of the interaction of non-stationary shock waveswith barriers. Z. Vycisl. Mat. i Mat. Fiz., 1:267–279, 1961.

[35] C. Zhang and R. J. LeVeque. The Immersed Interface Method for acoustic waveequations with discontinuous coefficients. Wave Motion, 25:237–263, 1997.


Date post:	10-Apr-2020
Category:	Documents
Upload:	others
View:	0 times
Download:	0 times

The coupling of homogeneous models for two-phase...

Documents