STUDY OF RELAXATION PROCESSES IN A TWO-PHASE FLOW MODEL

. This work concerns the analysis of the relaxation processes toward thermodynamical equilibrium arising in a compressible immiscible two-phase flow. Classically the relaxation processes are taken into account through dynamical systems which are coupled to the dynamics of the flow. The present paper compares two types of source terms which are commonly used: a BGK-like system and a mixture entropy gradient type. For both systems, main properties are investigated (agreement with second principle of thermodynamics, existence of solutions, maximum principle,...) and numerical experiments illustrate their asymptotic behaviour.


Introduction
This paper enters the framework of the modelling of compressible multiphase flows. Numerous models have been developed in order to depict their dynamics as well as the thermodynamical disequilibrium between the phases [1,2,16]. several approaches exist, we focus here on the two-phase modelling. On the one hand, the two-fluid models consider that the two phases evolve with their own velocities, leading to a model composed of two Euler-type systems coupled through nonconservative terms involving so-called interfacial quantities and relaxation source terms [1]. In this paper, we concentrate on the one-velocity approach, so that the disequilibrium only concerns thermal and mechanical exchanges as well as mass transfer between the two phases. Such models, referred as homogeneous models [2,10,15,18], enjoy many pros (compared to two-fluid models), notably a conservative form and thus well defined schock relations. In the homogeneous approach, one can assume that the phases are at thermodynamical equilibrium, leading to the class of Homogeneous Equilibrium Models. In this case, the system corresponds to one Euler system (for the mixture conserved quantities) with a complex closure pressure lawP (τ, e) which has to depict the thermodynamical behaviour of the mixture, for a given state of volume τ and internal energy e [4]. Such models are often inoperable because of the complexity of the mixture pressure law. A convenient way to get rid of this complexity is to relax the pressure law by making it depend on additional quantities Y , for instance the fractions of volume α, of mass y and energy z of one of the two phases. Then the pressure law is P (τ, e, Y ) and the system is completed with additional equations on these fractions with source terms Γ(τ, e, Y ), which have to depict the relaxation towards the thermodynamical equilibrium Y eq , such that P (τ, e, Y eq ) =P (τ, e). by the dynamical system dY dt (t) = Γ(τ, e, Y ).
Admissible source terms have to fulfill two major constraints: (1) Second principle of thermodynamics: the source term has to guarantee the growth of the thermodynamical entropy s(τ, e, Y ) of the mixture (and, similarly, the decrease of the mathematical entropy). A sufficient condition in the case of the Homogeneous Relaxation Models is to impose Γ(τ, e, Y ) · ∇ Y s(τ, e, Y ) ≥ 0; (2) Asymptotic states: for a given state (τ, e) of the fluid, the source terms have to capture the right thermodynamical equilibrium, that is The aim of the present paper is to compare two possible forms of source terms. The first one, recently referred as BGK source term [12], corresponds to a linearization around the equilibrium state Y eq . Its apparition goes back to [2,7,14] and is now used in multiphase configurations with complex phasic equation of state (EOS) [11,19]. The alternative source term corresponds to the gradient form of the mixture entropy, which directly fulfills (2) [18].
In Section 2 the thermodynamics of the two-phase mixture is presented, considering that both phases satisfy a Stiffened Gas law and are immiscible. The thermodynamical equilibrium is defined as the result of the maximization of the mixture entropy and its standard properties are recalled. Section 3 and 4 gather properties of these two choices of source terms, when considering a mixture of two Stiffened Gases. The advantage of this particular equation of state is that it allows to obtain relevant results while keeping a reasonable degree of complexity. In both sections, we investigate the eligibility and consistency of the source terms and list their major features: asymptotic state, trajectories, numerical resolution... Some remarks will be given on the influence of relaxation time scales. In the concluding Section 5, a comparison is proposed.

Thermodynamics of an immiscible two-phase mixture
This section concerns the thermodynamical description of the two-phase mixture. We adopt an intensive description, which can be deduced from an extensive one, see [5,15] for a proper derivation. We consider a fluid of intensive volume τ and internal energy e, composed of two phases that we label k = 1, 2. Each phase k is depicted by its phasic intensive volume τ k and internal energy e k . Following standard thermodynamics, the phase k is entirely described by an entropy function, which properties are stated in the following paragraph.
Then while investigating the mixture intensive entropy, we characterize the thermodynamical equilibrium. It results from a maximization of the mixture entropy under the constraints of mass conservation, energy conservation and immiscibility of the two phases. These three constraints are expressed in terms of fractions of mass y k , energy z k and volume α k of the phase k, so that the maximization process provides a characterization of the equilibrium vector Y eq . Note that the miscible case is treated in Appendix A.

Phasic equation of state
We focus here on the intensive description of the phase k: a unit of mass can be described by its specific volume τ k > 0 and its specific energy e k > 0, using an entropy function (τ k , e k ) → s k (τ k , e k ). By adopting the Gibbs formalism, the entropy function complies with the following differential form where the phasic temperature T k and pressure P k are defined by The phasic chemical potential is defined by the relation The phasic entropy functions are supposed to be • strictly concave function of (τ k , e k ) on (R + * ) 2 , • of class C 2 , such that ∀(τ k , e k ) ∈ (R + * ) 2 , T k (τ k , e k ) > 0. These standard assumptions derive from hypothesis usually made on the extensive entropies, see [15].
In the sequel, we will focus on a particular choice of phasic equation of state, namely the Stiffened Gas equation. This equation has the advantage of being a quite realistic model while keeping a reasonable degree of complexity in the following analysis.
The phasic entropy function is defined by where C k > 0 is the heat capacity, −Π k is the minimal pressure, Q k is a reference enthalpy, γ k > 1 is the adiabatic coefficient and s 0 k > 0 is a reference specific entropy. This is an extension of the perfect gas law [17], which corresponds to the simpler case Π k = Q k = 0.
For further numerical experiments, we focus on parameters given in [19] for a liquid (with index 1) and vapor water (with index 2), which have been determined by an optimization process

Mixture entropy function
Let w = (τ, e) be the state vector of the two-phase fluid. It is composed of two phases of specific volume τ k and e k , such that where y k ∈ [0, 1] is the mass fraction of the phase k, α k ∈ [0, 1] the volume fraction and z k ∈ [0, 1] the energy fraction. The mass and energy conservation impose that while assuming that the two phases are immiscible, the volume constraint is Note that Appendix A investigates the case of a miscible mixture which corresponds to the constraint α 1 = α 2 = 1.
In order to simplify notations, we note Y = (y, α, z) = (y 1 , α 1 , z 1 ) such that The intensive entropy of the mixture is given by It is defined on a convex subset of (R + ) 2 ×]0; 1[ 3 . Following [15], in order to capture single phase configuration, namely Y = (0, 0, 0) (pure vapor) or Y = (1, 1, 1) (pure liquid), we extend the entropy function on with notations 0 R 3 = (0, 0, 0) and 1 R 3 = (1, 1, 1), by The expression (11) can be deduced from the extensive form of the mixture entropy, see [4,5,7,15]. As a consequence, concavity properties result from those of the extensive entropy, namely: • for a given Y , (τ, e) → s(τ, e, Y ) is concave, • for a given (τ, e), Y → s(τ, e, Y ) is concave. In the case of a mixture of two Stiffened Gases, the entropy s is Definition 2.1. The set of admissible fractions for a given intensive state w = (τ, e) ∈ (R + * ) 2 is Proposition 2.2. For any couple w = (τ, e) ∈ (R + * ) 2 , C w is an open unit cube subset bounded by two planes P 1 and P 2 : Figure 1 represents the set of admissible fractions C w for a given state (τ = 0.01, e = 10 7 ) and Stiffened Gas law s k , k = 1, 2, with parameters (1). The planes P 1 and P 2 are represented in red and blue. Observe that the pure states 0 R 3 and 1 R 3 belong to the planes P 1 and P 2 respectively. In the case of perfect gases, the plane P 1 (resp. P 2 ) coincides with the plane {z = 0} (resp. {z = 1}).  (1). The state 0 R 3 belongs to the plane P 1 (in red) and the state 1 R 3 to the plane P 2 (in blue), see Proposition 2.2. In this particular case, one remarks that E ∈ C w while B, C, G, I, J ̸ ∈ C w .

Thermodynamical equilibrium
The second law of thermodynamics states that for any thermodynamical evolution, the entropy of the mixture cannot decrease and that the equilibrium corresponds to the maximum of the entropy. Definition 2.3. For a given w = (τ, e) ∈ (R + * ) 2 , the equilibrium fraction Y eq ∈ C w maximises the mixture entropy Y → s(τ, e, Y ) on C w . In other words The study of this optimization procedure has been the topic of numerous studies, notably for immiscible mixtures of Stiffened Gas law [4,10,11,19].
When the maximum is reached in the interior of the domain C w , then the equilibrium corresponds to a saturation state [7,14].
Proposition 2.4. When the maximum of the entropy is reached in the interior of the set of constraints C w , it is characterized by the equality of phasic pressures, temperatures and chemical potentials: Let us point out major features.
• When the maximum is reached in the interior of C w , it is not sure that Y eq is well defined. Indeed, the entropy function is only concave (and not strictly concave) on C w , so that Y eq is not necessarily uniquely defined, see Remark 2.5.
• The continuous extension of the mixture entropy guarantees that Y eq can take the values 0 R 3 or 1 R 3 , that correspond to the single phase equilibria. • The difficulty of the computation of Y eq directly depends of the chosen equations of state. For a mixture of two Stiffened Gases, it consists of a algebraic problem, given in [10], but for more complex EOS such as Noble-Able-Stiffened Gases (NASG) or NASG-Chemkin laws [19], non explicit changes of variables lead to tedious computations, which become even worse when considering additional phases, see [13,19]. • In [9], the authors affirm that the equilibrium is a single phase state in most cases. It would be interesting to verify this on a statistical study by choosing a relevant domain for (τ, e).
Remark 2.5. Attempts of proof concerning this strict concavity of the mixture entropy have been provided for instance in [13]. They rely on the study of the concavity of the extensive mixture.
The idea is to show, for a given extensive state W = (M, V, E) (with mass M , volume V and energy E), that the extensive mixture entropy S is a strictly concave function when restricted on the set H(M ) of the states whose total mass is equal to M . If this property is true, it follows the strict concavity of the intensive entropy.
In order to prove the strict concavity of S on H(M ), the idea is to show that the degeneracy manifold of S, ker ∇ 2 W S(W ) coincide with Vect(W ), that is to say: If so, then the restriction of [15]. One remarks that ker ∇ 2 W S(W ) is a vectorial subset in the 6-dimensional space in the case of a two-phase flow. The second inclusion seems difficult to show.
Since the strict concavity property of the mixture entropy s is crucial to model a physically relevant system, a brief numerical study has been done on the sign of ∇ Y s(τ, e, Y ). No fraction such as ∇ Y s(τ, e, Y ) = 0 have been found for the case presented on Figure 1, for 10 9 values tested in ]0; 1[ 3 .

BGK-like source term
This first type of source terms corresponds to a BGK linearized form for two-phase mixture. It has been introduced in the homogeneous literature, first in the two-phase case [2,7,10,14] and more recently in the multiphase context [13,19].
For a given state (τ, e, Y ) ∈ C, supposing Y → s(τ, e, Y ) is strictly concave, there exists a unique equilibrium fraction state Y eq which maximises the entropy s, where Y eq only depends of (τ, e), see Definition 2.3.
This definition involves time scales λ which generally depend on time and are distinct. They correspond respectively to the mass transfer, the mechanical and thermal relaxation times.
Proof. It is obvious that the stationary states of the source term Γ 1 are vector Y eq of Definition 2.3. It remains to check the entropy growth. Let (τ, e, Y ) ∈ C be a given state. It holds thanks to the concavity of Y → s(τ, e, Y ). □ Note that invoking the concavity argument requires to consider a unique relaxation parameter λ. Recently a modification of BGK-like source terms which allows distinct time scales has been introduced in [12].
admits a global solution Y given by where Λ is a primitive of 1/λ. In the case of a constant in time relaxation parameter λ, it holds This last expression can be interpreted as a barycenter of the initial state Y (0) and the equilibrium state Y eq . Trajectories of the solutions are exponential towards equilibrium fractions. Figure 2 illustrates this behaviour for a mixture of two Stiffened Gases with parameters (1). One observes the evolution in time of the fractions Y (t) for an initial state Y (0) = (0.5, 0.5, 0.5), with specific volume τ = 5.5.10 −3 and internal energy e = 10 8 . The time scale is constant in time and set to λ = 1. The equilibrium fractions, represented by the dotted lines, are reached at time t = 6s. Another important point is that the time when the asymptotic state is reached does not depend on the initial state of the fluid.
One major advantage of the BGK-like source term is that it can be explicitly integrated. However it requires to determine the equilibrium state Y eq for every considered state (τ, e, Y ) ∈ C, that is to say solve the maximization problem (14). For two-phase mixture, and especially Stiffened Gas law, the maximization can be done explicitly [10]. But when more complex EoS are considered and the number of phases increases, the maximization has to be realized with an appropriate optimization algorithm: one is proposed in the case of three Stiffened Gases in [3], another one uses a Broyden method for a mixture of two stiffened gases plus a NASG-Chemkin law in [19]. When coupled with the fluid dynamics, for instance through a time splitting approach, this maximization procedure has to be performed in each control volume and at each time step, leading to a complex resolution and higher CPU costs, see [3].

Entropy gradient source term
The second type of source terms corresponds to the entropy gradient such that it directly satisfies the eligibility condition (2).
with λ = (λ 1 , λ 2 , λ 3 ) relaxation parameters. Using the mixture entropy definition (11) and the phasic potentials (5)-(6), Γ 2 rewrites In the case of a mixture of two Stiffened Gases, the definition domain of Γ 2 is C w , represented in Figure 1. The source term Γ 2 corresponds to the entropy gradient times Y (1 − Y ) (understood as term by term product). This modification is mandatory for numerical purposes. Indeed the set of definition of the mixture entropy Y → s(Y ) is C w which is included in ]0; 1[ 3 . Thus a solution Y of the dynamical system Y ′ (t) = ∇ Y s(τ, e, Y ) is such that Y (t) ∈]0; 1[ 3 , for all t > 0 (unless the dynamical system is not defined). However in numerical applications, this maximum principle is not guaranteed, even for smaller time steps, as illustrated in Figure  3-left. It corresponds to a mixture of two stiffened gases with parameters (1) depicted by the dynamical system Y ′ (t) = ∇ Y s(τ, e, Y ) with initial state Y (0) = (0.5, 0.5, 0.5), τ = 8.10 −4 , e = 10 8 and λ = 1. One observes that the volume fraction becomes nonpositive at time t < 0.0025s. The correction term Y (1 − Y ) ensures the maximum principle as illustrated in Figure 3-right. Of course, this correction term modifies the trajectories, in particular the single phase equilibrium state seems to be reached at a longer time in the Figure 3 example. Every approximated solutions of the dynamical system (21) have been obtained with Adams-Bashforth plus backward differenciation formula types algorithms, using Python Scipy library.   Eligibility and Consistency. Assume that there exists λ 0 ∈ R * + such that ∀t ≥ 0, ∀k = 1, 2, 3, λ k (t) < λ 0 . For any given state (τ, e) ∈ (R + * ) 2 , the source term Γ 2 complies with the entropy growth criterion (2).
which concludes the proof. □ Unlike the previous BGK-like source terms, solving exactly the system is out of reach. However existence en uniqueness of local solution to Cauchy problem is guaranteed by the Cauchy-Lipschitz theorem. The following purpose is to ensure that the stationary states of the dynamical system (1)-(21), namely vector Y such that Γ 2 (Ȳ ) = 0, correspond to physically relevant thermodynamical equilibria of the two-phase mixture, in the sense of (3).
On the other hand, since Γ 2 is not defined on the border of the cube [0, 1] 3 , one cannot determine if single phase statesȲ = 0 R 3 orȲ = 1 R 3 are stationary states of the dynamical system and if these states are attractive (using standard arguments like the jacobian linearization). In a similar way, (1, 0, 0) could be a stationary state of the dynamical system, even if it has no physical sense. The following result states that trajectories Y cannot reach the border of the cube [0, 1] 3 in the case of a mixture of two Stiffened Gases with parameters (1). Proof. The proof relies on vector field argument. We study the limits of Γ 2 in the neighbourhood of the different considered sets, by supposing that it is defined in these areas, that depends on the thermodynamical parameters.
Since k 2 ̸ = 0, there is no possible equilibrium around the face {y = 0}. This case is illustrated on Figure 4 for several initial fractions Y (0) ∈]0, 1[ 3 (blue dots) with trajectories converging towards an asymptotic state (red triangle) corresponding to a saturation state. Observe that, since k 2 < 0, the curves are decreasing functions of α in the neighbourhood of the face {y = 0}.
The initial data Y (0) ∈]0, 1[ 3 is represented by a blue dot and the corresponding asymptotic state is represented in red.
• ∀y ∈]0; 1[, lim α→1, z→0 Four edges remain problematic because of undetermined limits. Finally, we can exclude four among the six corners: • lim y→1, α,z→0 For the two remaining corners (1, 1, 0) and (0, 0, 1), two components of Γ 2 do not admit a limit or the limit is undetermined. □ We emphasize that these arguments do not apply for the two corners (1, 1, 0) and (0, 0, 1) and the two faces {z = 0}and {z = 1}. Some improvements may be done while considering specific parameters Q k and π k . For instance, the realistic parameters (1) exclude some of these last sets as one can see on Figure 1, since Γ 2 is not defined on these areas. In these cases, we could investigate what happens on the planes, where the same phenomena can occur. The same kind of arguments hold since being in the neighbourhood of a plane induce limits towards ±∞ for some components of Γ 2 .
When considering a perfect gas mixture, the proof holds true and the corners (1, 1, 0) and (0, 0, 1) and the two faces {z = 0} and {z = 1} are not reachable as well.
In practical simulations, we did not find any case of trajectories with asymptotic states belonging to these sets. Some trajectories are plotted in Figure 5 in the case of an immiscible mixture of two Stiffened Gases with parameters (1), using standard Python semi-implicit algorithm that are sufficient for a two Stiffened Gas mixture. All the curves correspond to a solution of the dynamical system (1)-(21) with the initial data Y (0) = (0.5, 0.5, 0.5) with λ = 1. The internal energy is e = 10 8 and different specific volumes are considered. The green curve correspond to τ = 1.10 −3 and asymptotically converges towards the state 0 R 3 , that is to the presence of the phase 2 only. The blue curve corresponds to the case τ = 1.10 −2 and converges towards the single phase asymptotic state 1 R 3 (only the phase 1 is present). The red curve is obtained with τ = 5.5.10 −3 and corresponds to a saturation state Y eq ∈]0, 1[ 3 satisfying (2.4). Figure 5. Examples of trajectories of the dynamical system Γ 2 with initial data Y (0) = (0.5, 0.5, 0.5) and e = 10 8 . The green curve is obtained with τ = 1.10 −3 and converges towards the single phase state 0 R 3 . The blue curve corresponds to τ = 1.10 −2 and converges towards the single phase state 1 R 3 . The red curve with τ = 5.5.10 −3 converges to a saturation state in the interior of the cube. Figure 6 illustrates possible asymptotic states of the dynamical system (1)-(21). For a given internal energy e = 10 8 , each trajectory corresponds to a given initial state (τ, e, Y (0)), where Y (0) and τ are determined randomly respectively in ]0, 1[ 3 and in ]5.10 −4 ; 1.10 −2 [. The initial fraction Y (0) is represented by a blue dot and the asymptotic states with final time t = 0.1 are plotted in red. One observes that the asymptotic states are distributed in the interior of the cube and the corners 0 R 3 and 1 R 3 appear to be accumulation points.

Concluding comparison of the source terms
The BGK-like source terms provide very simple solutions, as soon as the equilibrium state Y eq is known for a given state (τ, e) ∈ (R + * ) 2 . The trajectories are of exponential type. In particular the time to return to equilibrium is always the same, no matter is the initial data. As the source term Γ 2 is concerned, we cannot give an explicit form of the solutions. In the general case, its trajectories are far more complex, notably non monotone, and the return to equilibrium directly depends on the initial state. Figure 7 presents trajectories for a two stiffened gas mixture with parameters (1) . They are obtained, for both source terms, for a state τ = 5.5.10 −3 , e = 10 8 with a relaxation time λ = 1. The same initial state Y (0) = (0.5, 0.5, 0.5) is considered. One observes that the trajectories for the BGK-like system have the expected exponential behaviour. On the contrary, trajectories of the entropy gradient source term are non monotone. An interesting point is the time when the asymptotic state is reached: the source term Γ 2 achieves the thermodynamical equilibrium state faster than the BGK-like source term.
To finish a comparison in terms of numerical complexity is mandatory. The BGK-like source term requires to compute the equilibrium fractions Y eq for any state (τ, e). When complex EoS are used, or if we consider more phases, optimization algorithm have to be implemented, see [3,11,19]. When coupled with the fluid dynamics, such algorithms are very CPU consuming since they are called in every cell at every time step.
Such complexity could be avoided using the source term Γ 2 . Of course, when considering complex phasic EoS or a larger number of phases, such source term can be stiff but semi-implicit or high-order explicit time integration method can be used and provide expected asymptotic states.
Finally, the entropy gradient source term allows to easily consider with different relaxation times for the mechanical, thermal exchange or the mass transfer. This is not the case for the BGK-like source term which requires the maximization of the mixture entropy with a common relaxation time, even if some recent improvements have been done on this topic [12].  Figure 7. Comparison of the trajectories of models Γ 1 and Γ 2 for the same initial data Y (0) = (0.5, 0.5, 0.5), τ = 5.5.10 −3 , e = 10 8 , λ = 1. While Γ 1 's trajectories are increasing exponential curves, Γ 2 's are not monotone. For the source term Γ 2 , the asymptotic state is reached at time t = 0.05s, and approximatively at time t = 3.5s for the source term Γ 1 . Asymptotic states are the same.
function (ρ k , ϵ k ) →s k (ρ k , ϵ k ). Note thats k = ρ k s k . We have the following relations, see [5] for a more complete description with the definitions and the pressure verifies Similarly to the specific case, we suppose that the phasic volumic entropies are strictly concave functions of their variables.
In the miscible case, both phases share the same volume, which gives equal volume fractions α 1 = α 2 = 1 Consequently, we only have two constraints, namely mass and energy conservation Hence the description in the miscible case requires one variable less than in the immiscible one. A unit of volume is described byw = (ρ, ϵ), where ρ is the volumic mass and ϵ is the volumic energy. We note in the following Y = (y, z) = (y 1 , z 1 ) the fraction of mass and energy of the phase 1. Finally, the volumic mixture entropys is given bỹ whose definition domainC is a convex subset of (R + ) 2 × (]0; 1[ 2 ∪0 R 2 ∪ 1 R 2 ), see [5] for a derivation from the extensive setting. Sinces is the sum and a composition of strictly concave and affine functions, we have the following results • for a givenỸ , (τ, e) →s(τ, e,Ỹ ) is strictly concave, • for a given (τ, e),Ỹ →s(τ, e,Ỹ ) is strictly concave.
When considering a miscible mixture of two Stiffened Gases, the definition domainCw ofỸ →s(ρ, ϵ,Ỹ ) is a subset of ]0; 1[ 2 bounded by two lines Following the Section 2.3, we can characterize the thermodynamical equilibrium of a miscible two-phase mixture. We emphasize two main differences: firstly, the entropy maximization is done in the two-dimensional space, and secondly, the equilibrium fractions are well defined thanks to the strict concavity ofỸ →s(τ, e,Ỹ ).
The entropy gradient source terms can be defined in the same way that in the immiscible case, but once again in the two-dimensional space. It reads with λ = (λ 1 , λ 2 ) relaxation parameters. Using the mixture entropy definition (26) and the phasic potentials (23)-(24),Γ 2 rewritesΓ .
Let us give similar properties of Proposition 4.4. It consists of studying the trajectories around the border of the square ]0; 1[ 2 . As in the immiscible case, the proof relies on vector field argument. The asymptotic behaviour ofΓ 2 is investigated in the neighbourhood of the edges of the square.
One remarks that this property can be stated for both edges {z = 0} and {z = 1} by choosing thermodynamical parameters and using the constraints (27).