Journal of Computational Physics

27  Download (0)

Full text

(1)

Contents lists available atScienceDirect

Journal of Computational Physics

www.elsevier.com/locate/jcp

A mixture-energy-consistent six-equation two-phase numerical model for fluids with interfaces,

cavitation and evaporation waves

Marica Pelanti

a,

, Keh-Ming Shyue

b

aDepartment of Mechanical Engineering, École Nationale Supérieure de Techniques Avancées – ENSTA ParisTech, 828, Boulevard des Maréchaux, 91762 Palaiseau Cedex, France

bDepartment of Mathematics, National Taiwan University, Taipei 106, Taiwan

a r t i c l e i n f o a b s t r a c t

Article history:

Received 31 March 2013

Received in revised form 30 October 2013 Accepted 1 December 2013

Available online 10 December 2013

Keywords:

Multiphase compressible flow models Mechanical relaxation

Thermo-chemical relaxation Cavitation

Phase transition Finite volume schemes Wave propagation algorithms Riemann solvers

We model liquid–gas flows with cavitation by a variant of the six-equation single-velocity two-phase model with stiff mechanical relaxation of Saurel–Petitpas–Berry (Saurel et al., 2009)[9]. In our approach we employ phasic total energy equations instead of the phasic internal energy equations of the classical six-equation system. This alternative formulation allows us to easily design a simple numerical method that ensures consistency with mixture total energy conservation at the discrete level and agreement of the relaxed pressure at equilibrium with the correct mixture equation of state. Temperature and Gibbs free energy exchange terms are included in the equations as relaxation terms to model heat and mass transfer and hence liquid–vapor transition. The algorithm uses a high-resolution wave propagation method for the numerical approximation of the homogeneous hyperbolic portion of the model. In two dimensions a fully-discretized scheme based on a hybrid HLLC/Roe Riemann solver is employed. Thermo-chemical terms are handled numerically via a stiff relaxation solver that forces thermodynamic equilibrium at liquid–vapor interfaces under metastable conditions. We present numerical results of sample tests in one and two space dimensions that show the ability of the proposed model to describe cavitation mechanisms and evaporation wave dynamics.

©2013 Elsevier Inc. All rights reserved.

1. Introduction

The modeling of cavitating flows is relevant in numerous areas of engineering, from naval and submarine systems design to aerospace and nuclear power plants technologies. Cavitating fluids are multiphase mixtures that often involve complex hydrodynamic and thermodynamic processes: liquid–vapor phase transition, dynamical creation of interfaces, vapor struc- tures collapse, and associated shock wave formation and interaction (cf.[1–3]). As a further reason of complexity, in many industrial applications these flows occur in irregular geometries and they have a multi-dimensional character.

Extensive work has been dedicated in the past decades to the simulation of cavitating flows and liquid–vapor flows with phase change, see for instance[4–15]and the references therein. Among the different modeling approaches, the class of hyperbolic compressible multiphase models stemming from the original model of Baer–Nunziato[16] has shown great capabilities in describing the complex wave patterns and thermodynamic mechanisms of cavitation. A first essential feature

*

Corresponding author. Tel.: +33 1 69 31 98 19; fax: +33 1 69 31 99 97.

E-mail addresses:marica.pelanti@ensta-paristech.fr(M. Pelanti),shyue@ntu.edu.tw(K.-M. Shyue).

0021-9991/$ – see front matter ©2013 Elsevier Inc. All rights reserved.

http://dx.doi.org/10.1016/j.jcp.2013.12.003

(2)

of these models is that compressibility is taken into account for all phases, vapor as well as liquid. This is fundamental to correctly capture wave propagation phenomena and acoustic perturbations, and it is particularly crucial when liquid–

vapor transition occurs[8]. Another important property is that these models can retain temperature and Gibbs free energy non-equilibrium effects, thus they are able to capture metastable states as well as evaporation fronts, when heat and mass transfer processes are included in the physical description through thermal and chemical relaxation source terms.

There exist various formulations of compressible temperature non-equilibrium multiphase flow models, depending on the assumptions on mechanical and kinetic phase equilibrium. In choosing a particular model, one has to find a good compromise between the accuracy of the description of the physical phenomena and the ability of conceiving robust and efficient numerical methods. In the present work, we are interested in the hyperbolic single-velocity six-equation model proposed by Saurel, Petitpas and Berry in[9]for compressible two-phase flows, see also Zein et al.[17]. This model consists of an advection equation for the volume fraction of one phase, mass and internal energy equations for each phase, and a mixture momentum equation. The six-equation model assumes instantaneous velocity equilibrium between the two phases, but it retains mechanical, thermal and chemical non-equilibrium effects. In the limit of instantaneous pressure relaxation the model reduces to the well known compressible two-phase flow model of Kapila et al.[18]. Nonetheless, as emphasized in[9], and as we briefly recall in Section2, numerically it is more advantageous to solve the six-equation system with stiff mechanical relaxation rather than the Kapila et al.[18]pressure-equilibrium five-equation model system.

The single-velocity six-equation two-phase model with stiff pressure relaxation was employed in [9]for applications to interface problems and mechanical cavitation processes (that is cavitation with no phase transition). It was later used by Zein et al. in [17] to simulate liquid–vapor transition in metastable liquids. One difficulty of the numerical algorithm illustrated in the latter work, as noted by the authors, is that it may require a very small time step for stability for some expansion problems with phase transition, due to the stiffness of the chemical relaxation terms. Only one-dimensional numerical results are presented by the authors in[17].

The aim of the present paper is to conceive a new multiphase flow computational model on the basis of the six-equation system of[9]that could deal efficiently with interfaces, cavitation and evaporation waves, while retaining simplicity and time-affordability. The key idea of our approach is to employ an alternative mathematical formulation of the standard six-equation model system [9]in the numerical discretization. Rather than using the two phasic internal energy equations of the classical model, in our algorithm we employ two equations for the phasic total energies. Mathematically, these two model systems are equivalent. The present model, however, is numerically advantageous with respect to the standard one, since it allows us to easily design a simple numerical method that ensures important consistency properties with mixture total energy conservation and with the mixture thermodynamic state. More specifically, first, we are able to automatically recover a conservative discrete form of the mixture total energy equation, whereas the classical six-equation model system needs to be augmented with an additional conservation law for the mixture total energy to correct the thermodynamic state [9,17]. Secondly, as a consequence of the mixture total energy conservation consistency property, we are able to easily ensure agreement of the relaxed pressure at equilibrium with the correct mixture equation of state for the full six-equation two-phase model that includes mechanical and thermo-chemical stiff relaxation effects. Relaxation terms are therefore efficiently handled.

To numerically solve the proposed two-phase model with pressure, temperature, and Gibbs free energy relaxation, we employ a simple fractional step approach that consists of the homogeneous hyperbolic system solution step, and a sequence of steps thereafter to solve systems of ordinary differential equations containing the relaxation source terms. A high- resolution wave propagation method based on Riemann solvers (HLLC and Roe) (cf. [19]) is employed for the numerical solution of the homogeneous hyperbolic system. The algorithm is easily implemented in the framework of the CLAWPACK software package[20]. For solving the ordinary differential equations with stiff relaxation sources, we have devised robust solvers that drive the mixture to the desired equilibrium conditions in a sequence of relaxation processes (cf. [21,10,8,17, 14]). In this procedure, similar to[8,22], thermodynamic equilibrium is forced at liquid–vapor interfaces under metastable conditions. Numerically for this task we employ an idea similar to[14,23]that uses the thermodynamic equilibrium condi- tions to reduce the solution of the ODEs relaxation problem to the solution of a simple system of algebraic equations for the equilibrium state variables.

This paper is organized as follows. In Section2.1, we begin by recalling the six-equation single-velocity model with stiff mechanical relaxation of Saurel–Petitpas–Berry[9]for compressible two-phase flows. We then propose in Section2.2a vari- ant of this model system, by employing phasic total energy equations in the mathematical formulation instead of the phasic internal energy equations of the classical approach. The extended model that includes thermal and chemical relaxation terms to model heat and mass transfer is described in Section2.3. In Section3we illustrate the numerical method to solve the basic model system with mechanical relaxation only. In this section we also discuss the mixture-energy-consistency property of the algorithm. The numerical treatment of temperature and Gibbs free energy relaxation source terms is de- scribed in Section 4. Finally, in Section5we present a selection of numerical results obtained by employing the proposed method with and without activation of heat and mass transfer.

(3)

2. The six-equation single-velocity two-phase flow model 2.1. Phasic-internal-energy-based formulation

The six-equation single-velocity compressible two-phase flow model with stiff mechanical relaxation proposed by Saurel et al.[9]has the form

t

α

1

+ 

u

· ∇ α

1

= μ (

p1

p2

),

(1a)

t

( α

1

ρ

1

) + ∇ · ( α

1

ρ

1u

 ) =

0

,

(1b)

t

( α

2

ρ

2

) + ∇ · ( α

2

ρ

2u

 ) =

0

,

(1c)

t

( ρ 

u

) + ∇ · ( ρ

u

 ⊗ 

u

) + ∇( α

1p1

+ α

2p2

) =

0

,

(1d)

t

( α

1E1

) + ∇ · ( α

1E1u

 ) + α

1p1

∇ · 

u

= − μ

pI

(

p1

p2

),

(1e)

t

( α

2E2

) + ∇ · ( α

2E2u

 ) + α

2p2

∇ · 

u

= μ

pI

(

p1

p2

).

(1f) Here

α

k is the volume fraction of phase k, k=1,2 (

α

1+

α

2=1),

ρ

k is the phasic density, pk is the phasic pressure, and Ek=

ρ

k

ε

k is the phasic internal energy, with

ε

k denoting the phasic specific internal energy. We have

ρ

=

α

1

ρ

1+

α

2

ρ

2 denoting the mixture density, andu representing the flow velocity vector. The source terms appearing in (1a),(1e), and(1f) model mechanical relaxation. In these terms

μ

>0 represents the pressure relaxation parameter and pI is the interface pressure, pI= Z2pZ11++ZZ12p2, where Zk=

ρ

kck2is the acoustic impedance of phase k, and ck is the sound speed of phase k. We assume an infinite-rate pressure relaxation with

μ

→ ∞, therefore mechanical equilibrium is reached instantaneously.

It is known (cf. [9]) that in this instantaneous pressure relaxation limit the six-equation model above reduces to the single-velocity single-pressure five-equation model of Kapila et al.[18](see also e.g.[24]). The five-equation model system is composed of two phasic mass balance equations, the mixture momentum equation, the mixture total energy equation, and an evolution equation for the volume fraction of one phase with a source term that results from the asymptotic limit of instantaneous velocity and pressure equilibrium of the non-equilibrium compressible two-phase flow model of Baer and Nunziato[16]. Despite this equivalence, the six-equation model(1)offers significant advantages in the numerical approxi- mation, compared to the five-equation model[18]. The main numerical issues in the solution of the five-equation system come from the non-conservative contribution in the volume fraction equation that depends on the divergence of the flow velocity and on the phasic impedances [9,59]. The variation of the volume fraction across acoustic waves associated to this term makes the construction of approximate Riemann solvers more challenging. In particular, the presence of this non-conservative contribution makes it difficult to preserve volume fraction positivity, especially when shocks and strong rarefaction waves are involved[9].

2.1.1. System’s closure

The closure of system(1)is obtained through the specification of an equation of state for each phase, which we choose to express in terms of Ek and

ρ

k, pk=pk(Ek,

ρ

k), k=1,2. The phasic sound speed can be written as ck=

κ

khk+

χ

k, where hk=Ekρ+kpk is the phasic specific enthalpy,

κ

k=pk(EEkk,ρk) and

χ

k=pk(Eρkk,ρk). The mixture sound speed of this model is

c

= 

Y1c21

+

Y2c22 (2)

where Yk=αkρρk is the mass fraction of phase k (Y1+Y2=1). Notice the monotonic character of(2)with respect to the volume fraction in contrast to the non-monotonic behavior of the Wood’s sound speed[25]of the Kapila et al. model[18].

This feature of the 6-equation model also represents an advantage over the 5-equation model in the numerical approxima- tion. As explained in [26], the non-monotonic behavior of the 5-equation model’s sound speed in the numerical diffusion zone of an interface may result in the presence of two sonic points in this region even when the flow is subsonic in the two pure fluids. The construction of approximate Riemann solvers able to handle robustly and efficiently these sonic points does not appear a simple task. In [9]the authors also explain that the sound speed non-monotonicity might affect the propagation of acoustic waves interacting with the interfacial zone, and result in a temporal delay in the wave transmission.

Let us remark that the pressure law pk=pk(Ek,

ρ

k), by using a more rigorous terminology, represents an incomplete equation of state. This law suffices to determine the fluid dynamics when thermal and chemical phenomena are neglected, as in the model above. However, a caloric law Tk=Tk(pk,

ρ

k)for each phasic temperature Tk is also needed in order to completely describe the thermodynamic state of the fluid (see e.g.[27]). The full thermodynamic characterization is required when heat and mass transfer terms are included in the flow model for problems with phase change, see Section2.3.

Here we will restrict our study to the case of species governed by the stiffened gas equation of state (SG EOS):

pk

(

Ek

, ρ

k

) = ( γ

k

1

)

Ek

γ

k

π

k

− ( γ

k

1

) η

k

ρ

k

,

(3a) Tk

(

pk

, ρ

k

) =

pk

+ π

k

Cvk

ρ

k

( γ

k

1

)

(3b)

(4)

for k=1,2, where

γ

k,

π

k,

η

k, and Cvkare material-dependent parameters. The associated entropy sk(pk,Tk)and the Gibbs free energy (chemical potential) gk(pk,Tk)=hkTkskare given by

sk

(

pk

,

Tk

) =

Cvklog Tkγk

(

pk

+ π

k

)

γk1

+ η

k

,

(3c)

gk

(

pk

,

Tk

) = 

γ

kCvk

η

k



Tk

CvkTklog Tkγk

(

pk

+ π

k

)

γk1

+ η

k

,

(3d)

where

η

k is a constant. With the SG EOS assumption(3), we have

κ

k= (

γ

k1),

χ

k= −(

γ

k1)

η

k, and the phasic sound speed can be expressed as ck=

γ

kpk+πk

ρk .

The mixture specific internal energy for the model considered here is defined as

ε

=Y1

ε

1+Y2

ε

2, and, equivalently, the mixture internal energy per unit volume is E =

ρε

=

α

1E1+

α

2E2. The latter relation, by using the isobaric assumption p1=p2=p in the energy laws Ek(pk,

ρ

k), k=1,2, gives the mixture equation of state, which determines implicitly the mixture pressure law p=p(E,

ρ

1,

ρ

2,

α

1):

E

= α

1E1

(

p

, ρ

1

) + α

2E2

(

p

, ρ

2

).

(4)

In the case with the SG EOS(3)we find an explicit expression for the mixture pressure p [28]:

p

(

E

, ρ

1

, ρ

2

, α

1

) =

E

− ( α

1

ρ

1

η

1

+ α

2

ρ

2

η

2

) − (

αγ11γ1π11

+

αγ22γ2π12

)

α1

γ11

+

γ2α21

.

(5)

When physically relevant values of the flow state variables are defined in the region of thermodynamic stability, the single- velocity six-equation model (1) is hyperbolic, that is it has real eigenvalues and a complete set of eigenvectors, see for instance [17]andAppendix A. Let us also recall that(1)represents the asymptotic limit of the hyperbolic seven-equation two-velocity two-phase flow model of Saurel and Abgrall [21] for instantaneous kinetic relaxation (see e.g. [17,18]). One advantage of the six-equation model with respect to the two-velocity model is that the order of the eigenvalues is a priori known, yielding an easy decomposition of waves in computing approximate solutions to Riemann problems (cf.[29–31]).

2.2. Phasic-total-energy-based formulation

For numerical reasons that we will discuss in the following, we propose to consider a mathematically equivalent for- mulation of the model (1), obtained by replacing the two phasic internal energy equations with two phasic total energy equations. We denote with Ek=Ek+12

ρ

ku· u the total energy of phase k. The alternative form of the six-equation model reads

t

α

1

+ 

u

· ∇ α

1

= μ (

p1

p2

),

(6a)

t

( α

1

ρ

1

) + ∇ · ( α

1

ρ

1



u

) =

0

,

(6b)

t

( α

2

ρ

2

) + ∇ · ( α

2

ρ

2



u

) =

0

,

(6c)

t

( ρ

u

 ) + ∇ · ( ρ

u

 ⊗ 

u

) + ∇( α

1p1

+ α

2p2

) =

0

,

(6d)

t

( α

1E1

) + ∇ · ( α

1E1u

 + α

1p1u

 ) + Σ(

q

,

q

) = − μ

pI

(

p1

p2

),

(6e)

t

( α

2E2

) + ∇ · ( α

2E2u

 + α

2p2u

 ) − Σ(

q

,∇

q

) = μ

pI

(

p1

p2

),

(6f) where the non-conservative termΣ appearing in the phasic total energy equations is

Σ (

q

,

q

) = −

u

· 

Y2

∇( α

1p1

)

Y1

∇( α

2p2

) 

= −

u

· 

(

Y2p1

+

Y1p2

)α

1

+ α

1Y2

p1

α

2Y1

p2

 ,

(6g)

with q denoting the vector of the system unknowns, see (6). Note that unlike the previous model (1) with the phasic internal energy equations, here the simple sum of the two non-conservative phasic total energy equations (6e) and(6f) recovers the equation expressing conservation of the mixture total energy E=E +12

ρ

u· u=

α

1E1+

α

2E2:

tE

+ ∇ · (

Eu

 + α

1p1u

 + α

2p2



u

) =

0

.

(7) This feature is beneficial in the numerical approximation of the model to ensure consistency with the conservation of E, see Section3.1.

(5)

2.3. Model with heat and mass transfer

A classical way to model thermal and chemical inter-phase phenomena mathematically consists in introducing additional heat and mass transfer source termsQandm, respectively, into the original model system (cf.˙ [14,32,21,8,11,17,22]). In the present work we follow in particular the modeling approach of Saurel and co-workers[8,22], employing the six-equation model with phasic total energy equations(6)as the basic system. The flow model augmented with heat and mass transfer terms takes the form:

t

α

1

+ 

u

· ∇ α

1

= μ (

p1

p2

) +

m

˙ ρ

I

,

(8a)

t

( α

1

ρ

1

) + ∇ · ( α

1

ρ

1u

 ) = ˙

m

,

(8b)

t

( α

2

ρ

2

) + ∇ · ( α

2

ρ

2u

 ) = − ˙

m

,

(8c)

t

( ρ 

u

) + ∇ · ( ρ

u

 ⊗ 

u

) + ∇( α

1p1

+ α

2p2

) =

0

,

(8d)

t

( α

1E1

) + ∇ · ( α

1E1



u

+ α

1p1u

 ) + Σ(

q

,

q

) = − μ

pI

(

p1

p2

) +

Q

+

eIm

˙ ,

(8e)

t

( α

2E2

) + ∇ · ( α

2E2



u

+ α

2p2u

 ) − Σ(

q

,

q

) = μ

pI

(

p1

p2

)

Q

eIm

˙ ,

(8f) where the termsQandm can be written as˙

Q

= θ(

T2

T1

),

(9a)

˙

m

= ν (

g2

g1

).

(9b)

Hereθand

ν

are the thermal and chemical relaxation parameters, respectively. These parameters are assumed to be infinite at selected locations, while they are set to zero elsewhere. More specifically, thermal transfer is activated at liquid–vapor interfaces, and thermo-chemical transfer is activated at liquid–vapor interfaces under metastable thermodynamic conditions (liquid temperature Tliqhigher than the saturation temperature Tsat at the given pressure):

θ =

 ∞

if



I

 α

1



1



I

,

0 otherwise

,

(10a)

ν =

 ∞

if



I

 α

1



1



Iand Tliq

>

Tsat

,

0 otherwise

.

(10b)

Here the parameter



Iidentifies liquid–vapor interface locations (e.g.



I=104). See[8]for further discussion on the choice of tolerances.

We refer the reader to[33,34]for a rigorous and systematic analytical study of Liu’s subcharacteristic condition[35]for the six-equation two-phase model with the three levels of relaxation considered here, namely, pressure relaxation, simultaneous pressure and temperature relaxation, simultaneous pressure, temperature and chemical potential relaxation. In particular, the works[33,34]show that the subcharacteristic condition (characteristic speeds of the relaxation system at least as large as the characteristic speeds of the relaxed equilibrium system) is satisfied for all cases as long as physical fundamental positivity conditions on a set of physically positive thermodynamic variables are fulfilled.

Let us notice that a theoretical pressure–temperature (p–T ) saturation curve can be obtained by imposing the equilib- rium condition on the Gibbs free energy g1=g2 for the liquid and vapor phases (cf.[36,37,8]). With the SG EOS(3a)–(3d), the p–T saturation curve is defined by the equation

AS

+

BS

T

+

CSlog T

+

DSlog

(

p

+ π

1

)

log

(

p

+ π

2

) =

0

,

(11a)

with

AS

=

Cp1

Cp2

+ η

2

η

1

Cp2

Cv2

,

BS

= η

1

η

2

Cp2

Cv2

,

CS

=

Cp2

Cp1

Cp2

Cv2

,

DS

=

Cp1

Cv1

Cp2

Cv2

.

(11b)

Here Cpk denotes the heat capacity at constant pressure, Cpk=Cvk

γ

k, k=1,2. The parameters in the SG EOS (3a)–(3d) can be chosen to fit the above theoretical p–T curve with the experimental curve for the considered material [36]. The theoretical p–T curve is used in the numerical algorithm, see Section4.2.

Concerning the interfacial density

ρ

Iand the interfacial specific total energy eIthat appear in the mass and heat transfer source terms, their expression can be obtained by imposing appropriate thermodynamic constraints on the thermal and chemical processes. For instance, in [17], Zein et al. determine these quantities in a model similar to (8) by assuming pressure and temperature equilibrium during the relaxation processes. However, thanks to the assumption (10) on the relaxation parameters and to the particular numerical algorithm used for the treatment of thermal and chemical source terms, there is no need to specify the expressions of

ρ

Iand eI(see Section4). Let us also note that, due to the symmetry of the source terms in(8), the mixture mass and total energy equations are recovered as for the model in(6).

(6)

To end this section, for the ease of later reference, we write(8)in a compact form as

tq

+ ∇ ·

f

(

q

) + σ (

q

,

q

) = ψ

μ

(

q

) + ψ

θ

(

q

) + ψ

ν

(

q

),

(12a) where

q

=

⎜ ⎜

⎜ ⎜

⎜ ⎜

⎜ ⎜

α

1

α

1

ρ

1

α

2

ρ

2

ρ

u

 α

1E1

α

2E2

⎟ ⎟

⎟ ⎟

⎟ ⎟

⎟ ⎟

,

f

(

q

) =

⎜ ⎜

⎜ ⎜

⎜ ⎜

⎜ ⎜

0

α

1

ρ

1u

 α

2

ρ

2u



ρ

u

 ⊗ 

u

+ ( α

1p1

+ α

2p2

)I α

1

(

E1

+

p1

)

u

α

2

(

E2

+

p2

)

u

⎟ ⎟

⎟ ⎟

⎟ ⎟

⎟ ⎟

, σ (

q

,

q

) =

⎜ ⎜

⎜ ⎜

⎜ ⎜

⎜ ⎜



u

· ∇ α

1

0 0 0

Σ (

q

,

q

)

−Σ(

q

,

q

)

⎟ ⎟

⎟ ⎟

⎟ ⎟

⎟ ⎟

,

(12b)

ψ

μ

(

q

) =

⎜ ⎜

⎜ ⎜

⎜ ⎜

⎜ ⎜

μ (

p1

p2

)

0 0 0

μ

pI

(

p1

p2

) μ

pI

(

p1

p2

)

⎟ ⎟

⎟ ⎟

⎟ ⎟

⎟ ⎟

, ψ

θ

(

q

) =

⎜ ⎜

⎜ ⎜

⎜ ⎜

⎜ ⎜

0 0 0 0

θ (

T2

T1

)

−θ(

T2

T1

)

⎟ ⎟

⎟ ⎟

⎟ ⎟

⎟ ⎟

, ψ

ν

(

q

) =

⎜ ⎜

⎜ ⎜

⎜ ⎜

⎜ ⎜

ν

g2ρg1

I

ν (

g2

g1

)

ν (

g2

g1

)

0

ν

eI

(

g2

g1

)

ν

eI

(

g2

g1

)

⎟ ⎟

⎟ ⎟

⎟ ⎟

⎟ ⎟

(12c)

withΣ as in(6g). Above we have put into evidence the conservative portion of the spatial derivative contributions in the system as∇ · f(q), and we have indicated the non-conservative term as

σ

(q,q). The source terms ψμ(q),ψθ(q), ψν(q) contain mechanical, thermal and chemical relaxation terms, respectively.

3. Numerical solution of the model system with mechanical relaxation

We begin by considering the solution of the model system(12)with only mechanical relaxation:

tq

+ ∇ ·

f

(

q

) + σ (

q

,

q

) = ψ

μ

(

q

).

(13) The treatment of thermal and chemical relaxation terms will be described in Section 4. To numerically solve the system above we use a fractional step technique, similar to [9,17], where we alternate between the solution of the homogeneous hyperbolic system and the solution of a system of ordinary differential equations that takes into account pressure relaxation source terms. That is, the algorithm consists of the following steps:

1. Homogeneous hyperbolic system. We solve over a time intervalt the homogeneous hyperbolic portion of(13):

tq

+ ∇ ·

f

(

q

) + σ (

q

,

q

) =

0

.

(14) 2. Stiff mechanical relaxation. We solve in the limit

μ

→ ∞the system of ordinary differential equations (ODEs)

tq

= ψ

μ

(

q

).

(15)

This step drives the two-phase flow to mechanical equilibrium with an equilibrium relaxed pressure p1=p2=p. In this step the partial densities, the mixture momentum, the mixture total energy, and the mixture internal energy remain constant. The volume fraction

α

1, the mixture pressure p, and the phasic internal energies

α

kEk,Ek=Ek(p, (

α

k

ρ

k)/

α

k) for k=1,2, are updated before returning to Step 1.

3.1. Mixture-energy-consistent discretization

Before illustrating each step of the algorithm in more details, we emphasize that in the design of a numerical method for computing approximate solutions of the two-phase model (12) it is important to ensure at all times conservation of the quantities that are physically conserved, namely the partial densities

α

k

ρ

k, k=1,2 (together with the mixture density

ρ

=

α

1

ρ

1+

α

2

ρ

2), the mixture momentum

ρ

u, and the mixture total energy E=

α

1E1+

α

2E2. Note in particular that the values of the equilibrium pressure to be used at the beginning of the homogeneous system solution step must satisfy the mixture EOS (4)for values of E =E12

ρ

u· u that correspond to conservation-consistent discrete values of E, so as to approximate correctly the flow thermodynamic state. Godunov-type schemes can be relatively easily formulated to preserve conservation at the discrete level of quantities that are governed by conservation laws (cf.[38,19]). However, the two-phase mathematical model (either in the form (1) or(6)) does not contain the conservation law for the total energy, but two phasic energy equations from which the total energy is recovered. The difficulty then is to discretize the phasic energy

(7)

equations in such a way that total energy conservation is fulfilled at the discrete level, and that consistency with the correct thermodynamic state is ensured in the sense made precise hereafter.

Let us denote with superscript 0 the quantities computed by solving the homogeneous system in Step 1 of the algorithm above, and with superscript∗the quantities at mechanical equilibrium computed in Step 2. As mentioned above (see also Section3.3), we have

ρ

k=

ρ

k0, k=1,2,(

ρ

u)= (

ρ

u)0,E=E0, and E=E0. Let us also denote with E0,C discrete values of the mixture total energy that come from a conservative approximation of the conservation law for E in(7).

Definition 3.1. We say that the numerical scheme based on the fractional step algorithm above is mixture-energy-consistent if the following two properties are satisfied

(i) Mixture total energy conservation consistency:

E0

=

E0,C

,

(16a)

where E0= (

α

1E1)0+ (

α

2E2)0. (ii) Relaxed pressure consistency:

E0,C

= α

1E1

p

, ( α

1

ρ

1

)

0

α

1

 + α

2E2

p

, ( α

2

ρ

2

)

0

α

2



,

(16b)

whereE0,C=E0,C(ρu)20ρ·(0ρu)0.

The first property (i) means that the sum of the discrete values of the phasic total energies given by the solution of the homogeneous system must recover discrete values of the mixture total energy that are consistent with a conservative discrete form of(7). The second property (ii) means that the value of the relaxed (equilibrium) pressure p predicted in the relaxation step must be equal to the pressure as computed through the mixture equation of state p(E0,C,

α

1,

ρ

1,

ρ

2), defined by(4). That is, with the SG EOS, this consistency condition reads (cf.(5))

p

=

E

0,C

− (( α

1

ρ

1

)

0

η

1

+ ( α

2

ρ

2

)

0

η

2

) − (

αγ11γ1π11

+

αγ22γ2π12

)

α1

γ11

+

γα221

.

(17)

The mathematical formulation of the two-phase model with the phasic total energy equations (6)easily allows us to satisfy both properties (i) and (ii). To ensure the property (i), it suffices in Step 1 to apply a standard conservative scheme to the conservative portion of the energy equations(6e)and(6f), that ist(

α

kEk)+ ∇ · (

α

kEku+

α

kpku), k=1,2, and to discretize symmetrically the non-conservative contributionΣ appearing there. In such a way, the sum of the discrete non- conservative energy equations recovers a conservative discrete form of the mixture energy equation, as a consequence of the cancellation of non-conservative discrete contributions. The fulfillment of mixture total energy conservation consistency then easily enables us to ensure also the property (ii), the agreement of the relaxed equilibrium pressure with the correct mixture equation of state. See the simple pressure relaxation procedure for Step 2 in Section3.3.

On the other hand, it appears difficult to obtain a mixture-energy-consistent scheme if we apply an analogous frac- tional step algorithm to the classical six-equation two-phase model (1). Although clearly both formulations (1) and (6) mathematically recover the conservation law for the mixture total energy, it seems hard to discretize the phasic internal energies equations(1e)and(1f)in a way that recovers a conservative discrete form of(7). Indeed, numerical models such as[9,17] built on the formulation (1)need to augment the six-equation system with the equation for E. The additional conservation law for E is solved through a standard conservative scheme to obtain consistent discrete values E0,C. These values are then used to correct the thermodynamic state predicted by the non-conservative internal energy equations via the mixture equation of state. Note that this approach in general does not guarantee the consistency property (ii), that is p=p(E0,C,

α

1,

ρ

1,

ρ

2).

3.2. Homogeneous hyperbolic system solution step

In Step 1 of the algorithm we employ the wave propagation method of[19,39]to compute approximate solutions of the homogeneous system(14). This method belongs to a class of Godunov-type finite volume schemes[40,41,38,19]for solving hyperbolic systems of partial differential equations. We describe hereafter the basic ideas of the method in one dimension.

The two-dimensional scheme will be briefly recalled in Section3.2.4.

3.2.1. One-dimensional wave-propagation scheme

We consider the solution of the one dimensional systemtq+ ∂xf(q)+

σ

(q, ∂xq)=0 (as obtained by settingu=u and

∇ = ∂x in(12)). We assume a uniform grid with cells of widthx, and we denote with Qin the approximate solution at the ith cell at time tn, i∈ Z, n∈ N. Settingt=tn+1tn, the second-order one-dimensional wave-propagation scheme has the form

(8)

Qin+1

=

Qin



t



x



A+



Qi1/2

+

A



Qi+1/2





t



x



Fih+1/2

Fhi1/2



,

(18)

whereAQi+1/2are the so-called fluctuations at interfaces xi+1/2between cells i and(i+1), and Fih+1/2are second-order correction fluxes for higher resolution.

Here the fluctuationsAQi+1/2are computed by solving local Riemann problems at cell interfaces xi+1/2 for each pair of data Qin, Qin+1corresponding to adjacent cells. A Riemann solver (cf.[41,38,19]) must be provided to perform this task.

Let us specify the jump relations that need to be satisfied by an approximate Riemann solver for a Riemann problem with left and right initial data qand qr. The Riemann solution structure defined by the solver can be expressed in general by a set of M waves Wl and corresponding speeds sl, M6. For example, see the HLLC-type solver (M=3) and the Roe-type solver (M=6) presented below. By using the formalism introduced in[42], we also define the f -wavesZl=slWl, l=1, . . . ,M, which have the dimension of a flux. The sum of the waves must be equal to the initial jump in the vector q of the system variables:



q

qr

q

=



M l=1

Wl

.

(19)

Moreover, for any variable of the model system governed by a conservative equation the initial jump in the associated flux function must be recovered by the sum of the associated f -wave components. In the considered model the conserved quantities are

α

k

ρ

k, k=1,2, and

ρ

u, therefore in order to guarantee conservation we need:



f(ξ )

f(ξ )

(

qr

)

f(ξ )

(

q

) =



M l=1

slWξl

=



M l=1

Zlξ (20)

for ξ=2,3,4, where f(ξ ) is the ξth component of the flux vector f , and Wlξ andZξl denote theξth component of the lth wave and of the lth f -wave, l=1, . . . ,M, respectively. It is clear that conservation of the partial densities ensures conservation of the mixture density

ρ

=

α

1

ρ

1+

α

2

ρ

2. In addition, we must ensure conservation of the mixture total energy, that is the consistency condition (i) in(16a):



fE

fE

(

qr

)

fE

(

q

) =



M l=1

sl



Wl5

+

W6l



=



M l=1

Zl5

+

Z6l

,

(21)

where fE=u(E+

α

1p1+

α

2p2)is the flux function associated to the mixture total energy E.

Once the Riemann solution structure {Wil+1/2,sli+1/2}l=1,...,M arising at each cell edge xi+1/2 is defined through a Rie- mann solver, the fluctuationsAQi+1/2 and the high-resolution correction fluxes Fih+1/2 in(18)are computed as

A±



Qi+1/2

=



M l=1



sl

i+1/2



±

Wil+1/2

,

(22)

s+=max(s,0), s=min(s,0), and Fhi+1/2

=

1

2



M l=1



sli+1/2



1



t



x



sli+1/2

 

Wih+,l1/2

,

(23)

whereWih+,l1/2are a modified version ofWli+1/2obtained by applying toWil+1/2a limiter function (cf.[19]).

We present in the next sections two approximate Riemann solvers that we have developed for the model system(12):

a HLLC-type solver[43,44,38], similar to the solvers described in[9,17], and a new Roe-type solver[45]. The conservation consistency condition(21)(equivalently(16a)) is easily fulfilled in both cases.

3.2.2. HLLC-type solver

To begin with, we define an approximate solver for the two-phase system(14)by applying the idea of the HLLC solver of Toro et al. [44] (see also[38]). One difficulty in designing a HLLC-type solver for the present model is related to the non-conservative character of the phasic energy equations, for which we lack a notion of weak solution in the distributional framework (see for instance on this subject[46–49]). Nonetheless, note that to correctly set the initial thermodynamic state for the solution of the homogeneous system we only require that the sum of the phasic total energies computed at the pre- vious time level fulfill the consistency condition(

α

1E1)0+ (

α

2E2)0=E0,C. The individual phasic energies are re-set at the beginning of Step 1 by using the relaxed variables p and

α

1, obtained in Step 2 through a procedure that by construction ensures that the energies’ sum recovers the correct mixture energy state (Section3.3). The individual phasic energies values that come from the solution of the homogeneous system have only a role in the initial condition for the mechanical relax- ation step. Aiming at designing the simplest method that could provide reasonable mixture-energy-consistent estimates for

Figure

Updating...

References

Related subjects :