VAPOUR-LIQUID PHASE TRANSITION AND METASTABILITY

. The paper deals with the modelling of the relaxation processes towards thermodynamical equilibrium in a liquid-vapour isothermal mixture. Focusing on the van der Waals equation of state, we construct an constrained optimization problem using Gibbs’ formalism and characterize all possible equilibria: coexistence states, pure phases and metastable states. Coupling with time evolution, we develop a dynamical system whose equilibria coincide with the minimizers of the optimization problem. Eventually we consider the coupling with hydrodynamics and use the dynamical system as a relaxation source terms in an Euler-type system. Numerical results illustrate the ability of the whole model to depict coexistence and metastable states as well.


Introduction
This paper deals with the modelling of vapour-liquid phase transition and the relaxation towards the thermodynamical equilibrium in liquid-vapour mixture. Focusing on isothermal problems, vapour bubbles may appear in a liquid through a decompression process, which leads to a pressure decrease below the saturation pressure and cavitation. If the process of depressurization is slow or weak enough, then the system may remain entirely liquid, although the pressure is below the saturation pressure. There is a delay in the phase transition, which may be shortened by a (strong enough) dynamic perturbation. This liquid state is commonly called metastable. A similar behavior holds in a vapour system submitted to a compression. All these phenomena result from the coupling between hydrodynamics and thermodynamics, the former inducing dynamical perturbations of thermodynamical equilibrium states.
Our purpose is therefore to provide a characterization of thermodynamical equilibria which on the one hand copes with coexistence (that is saturation of liquid and vapour), pure phases and metastable states, on the other hand can be coupled in a consistent way to the equations of hydrodynamics, in the present case the isothermal Euler system. In a first step, we obtain a static characterization of the equilibrium states of a fluid described by a single nonconvex energy, namely the van der Waals equation of state. Following Gibbs' formalism, this is described by a constrained minimization problem, introducing heterogeneity in the system. In doing so, we recover mathematical formulations of classical results in thermodynamics, such as the Gibbs phase rule and the Maxwell equal area construction. In a second step, we proceed towards the coupling to time evolution equations by providing an alternative characterization, precisely in terms of dynamical perturbations. Pure phases, metastable states and coexistence states are now attractive points of a carefully chosen dynamical system, with specific basins of attraction. In particular, metastability corresponds to some competition between two such basins.
In the first section of this paper we describe the Gibbs formalism for the van der Waals equation of state, and give the static characterizations of equilibrium states. The second section is devoted to the study of dynamical fluctuations around these states, with a particular focus on the precise description of the basins of attraction. Finally, in Section 3 we provide numerical illustrations of the coupling between thermodynamics and the Euler equations, following a strategy of relaxation with the thermodynamically consistent model previously developed.
The model we propose here is a variant of the one described in [8], and brings some precisions on the analysis of the dynamical system. We refer to this paper for complementary results, in particular concerning the coupling, and a more complete reference list.

Thermodynamical setting
This section is devoted to the description of the Gibbs formalism to describe thermodynamical equilibria. We show how to go from simple, homogeneous systems, which obviously cannot describe two-phase fluids, to more complex, heterogeneous systems. We investigate the notion of stable equilibria and give their characterization through the minimization of a nonconvex function.

Simple system modelling
Consider a single fluid with mass M > 0 and volume V > 0, assumed to be homogeneous and at rest, at a fixed temperature T . Then its description is completed by another quantity called the Helmholtz free energy F . According to the Gibbs' formalism, we say that the fluid is at equilibrium if its Helmholtz free energy is a function, also denoted by F , of its mass M and volume V : All the variables involved here are extensive, which means that they have the same scaling as the volume V : this corresponds to the notion of homogeneity of the sample. As a consequence, the function F has to be positively homogeneous of degree 1: Without loss of generality, we assume here that F is of class C 2 with respect to V and M . Two fundamental quantities can then be considered, namely the pressure p and the chemical potential µ, defined by Notice that these quantities are defined only when the system is at equilibrium. We recover the classical thermodynamic relation for isothermal fluids From homogeneity, we obtain the so-called Euler relation, which is known in this context as Gibbs' relation: An immediate computation shows that the pressure and the chemical potential are positively homogeneous of degree 0. In other words, they do not scale with the volume. Such quantities are called intensive quantities, typical other ones, which are defined out of equilibrium as well, are the specific volume τ = V /M , and the density ρ = M/V .
In the sequel, we shall work with intensive variables, so we rewrite all the preceding relations in terms of the density ρ. We first introduce the Helmholtz free energy per unit volume f Keeping the same notations p and µ, the pressure and the chemical potential as functions of the density ρ are given by Thus we obtain the intensive form of the Gibbs relation (5) together with the following relations A formula relating a thermodynamical quantity (e.g. the free Helmholtz energy or the pressure) to the mass and volume, or to the density, is called an equation of state (EoS). Usually equations of state are given in an incomplete way, for instance we only know the pressure law p(ρ), as in the usual perfect gases law: p(ρ) = RT ρ, T being the given temperature and R the gas constant. The formalism we use here requires the complete equation of state, that is the energetic representation (1). Recovering such a representation from a partial equation of state can be a challenge, see e.g. [1,2]. In the simple case of the perfect gas law, the function f is convex, which turns out not to be appropriate to represent phase transitions.
We shall use in the following the so-called reduced complete van der Waals equation of state. The van der Waals law is not a physically relevant EoS, but it enjoys the basic properties that allow a toy-model for phase transitions. It is the first example of so-called third order equations of state [13]. It consists of a family of functions of ρ, depending on the temperature T as a parameter. There exists a critical temperature T C above which this function is convex (the pressure is monotone), and below which it is not convex: there are two inflexion points at densities ρ − < ρ + , between which the pressure is a decreasing function. For T = T C , there exists a unique inflexion point, which defines a critical density ρ C and a critical pressure p C . Considering normalized critical quantities, that is setting ρ C = 1, p C = 1 and T C = 1, the reduced van der Waals EOS reads for 0 < ρ < 3. For more details on this classical nonconvex EoS, we can refer among many references to [1,9].

Inhomogeneity and equilibria
In the previous description we assumed homogeneity of the system. This is obviously no longer relevant to cope with phase transition, so that we introduce inhomogeneity into the model, and therefore have to investigate its thermodynamical stability. For given mass M and volume V , we split the system into an arbitrary number I ≥ 1 of simple subsystems, each one being described by the same EoS. Each subsystem has a mass M i , a volume V i and its free energy is given by F (M i , V i ), with the same F function, describing the fluid under consideration.
Assuming immiscibility of the subsystems, we recover total mass and volume by and the total free energy is given by . Among these decompositions, we have to select the physically admissible ones: this is the notion of thermodynamical stability. According to the Gibbs' formalism, the stable states are those which minimize the total Helmholtz free energy F under the constraints (11). Thus we have to solve It is more convenient to reformulate the problem in terms of intensive variables, introducing the partial densities A straightforward computation leads to the equivalent formulation of problem (12): Two remarks are in order here. First we do not take into account the positivity constraints on the densities and volume fractions here, we choose to treat them in a second step, see below. Next, and more important, the number of subsystems I itself is a varying parameter in the optimization problem. This last problem is actually related to the so-called Gibbs phase rule, which gives an upper bound to the number of possibly coexisting phases [1,10]. A mathematical analogue of this rule results from the analysis of the minimization problem, more precisely from Carathéodory's theorem which states that the set of minima for the problem (13) for a general function f in R n is a (possibly degenerate) simplex. In other words, for n = 1 as in our case, we have I ≤ 2, that is at most two subsystems may coexist, which are now called phases, namely the vapour phase and the liquid phase. For more details, see for instance [11] in the general context and [4,8] for the application to thermodynamics.
Following this, we consider now I = 2 and define the Helmoltz free energy of the system by Then the minimization problem rewrites Since the value of I is fixed now, we can recover the case I = 1 in (15) either by α i = 0 or ρ 1 = ρ 2 . Notice that if ρ 1 = ρ 2 , we can solve for α i in the constraints above: Without loss of generality, we can assume ρ 2 > ρ 1 , so that the positivity of α i amounts to ρ 1 ≤ ρ ≤ ρ 2 .
Introducing the Lagrange multipliers λ α and λ ρ corresponding to the constraints in (15), we obtain straightforwardly the optimality conditions for the minima in (15) For two states a and b we introduce the relative Helmholtz free energy of a with respect to b: Lemma 1.1. The equilibrium states are (1) Pure states: (2) Coexistence states: ρ 1 < ρ < ρ 2 , with α i given by (16) and which uniquely define the couple of saturation densities ρ * 1 < ρ − < ρ + < ρ * 2 , as well as the saturation pressure p * and saturation chemical potential µ * .
The above equality of pressures and chemical potentials (19) can be rewritten in two equivalent ways: • equality of the relative free energies: • Maxwell's area rule on the chemical potential, see [8]: The first property will be used to study the dynamical systems in the following section. The second one is another classical characterization of coexistence states. Indeed the van der Waals Helmoltz free energy is nonconvex in the so-called spinodal zone. In the pressure-density plane, for a temperature T below the critical temperature T C , the spinodal zone corresponds to the domain where the pressure decreases with the density. It is delimited by the densities ρ − and ρ + , see Figure 1-left. It is well-known that states belonging to the spinodal zone are unstable. A classical regularization of the van der Waals EOS consists in replacing the function f by its convex hull, or equivalently by its Fenchel biconjugate f * * , which turns out to be characterized precisely by the minimization problem (15) Once again, we refer to [11] for the general mathematical results and to [2,6] for applications in thermodynamics. The computation of f * * amounts to define the two saturation densities ρ * 1 < ρ * 2 satisfying (19), for the temperature T , and replace both pressure and chemical potential by the constants p * and µ * , for all states in between ρ * 1 and ρ * 2 , see Figure 1. These are the so-called Maxwell lines, and for the chemical potential the areas above and below this line are equal: this is the meaning of (21). These states correspond to coexistence states, and the set of all coexistence densities defined for all temperatures below T C is called the saturation dome. The metastable zones correspond to the states belonging to the saturation dome outside the spinodal zone. Observe that the metastable zones [ρ * 1 , ρ − ] and [ρ + , ρ * 2 ] correspond to monotone branches of the van der Waals pressure.
We emphasize that the necessary conditions in Lemma 1.1 include all equilibrium states, regardless of their stability. In particular, we recover in item (1) the complete van der Waals curve, including physically unstable states (spinodal zone), metastable states and pure stable states. To proceed further a classical way consists in studying the local convexity of F in (14). We prefer here, following [8], to introduce a relaxation towards the equilibrium states, by means of a dynamical system.

A revisited dynamical system and its attraction basins
To build the appropriate dynamical systems, we impose two basic criteria, namely • long-time equilibria coincide with the physically relevant equilibria in Lemma 1.1, • the free Helmholtz energy is dissipated along trajectories.
Such dynamical systems in some sense simulate the fluctuations in the thermodynamical system. Since the number of degrees of freedom in the constraints is 2, we need only to provide a dynamical system of dimension 2. Several choices are possible to define the set of extended variables. We propose here to use the partial densities (ρ 1 , ρ 2 ), which are naturally involved in the minimization problem. Even within this set, the choice of the dynamical system is not unique.

Defining the dynamical system
We set r = (ρ 1 , ρ 2 ), and want to use relation (16) to define α i (r) and thus reduce the number of variables to 2. As mentioned before, to ensure in addition the positivity of the volume fractions, we must assume in this section that 0 <ρ 1 ≤ ρ ≤ ρ 2 < 3 and ρ 1 < ρ 2 .
In this context, the Helmholtz free energy of the system becomes a function of r, still denoted by F: Using we obtain where f (ρ 2 |ρ 1 ) is the relative energy of ρ 2 with respect to ρ 1 , as defined in (18). Both F and ∇F are defined only for ρ 1 < ρ 2 , however we have the following result.
We emphasize here that the equilibrium statesr = (ρ, ρ) are valid for all values of ρ and go over the whole van der Waals curve. This is in accordance with the definition of thermodynamical equilibrium. To go further and identify the physically admissible equilibrium states, we must investigate their stability and attractivity. Proposition 2.3 (Attractivity). The equilibrium states are classified as follows: is an attractive point, •r i , i = 1, 2 are unstable hyperbolic points.
Finally, the case (2b) leads to .
To complete the study of equilibrium states, in particular to cope with the degenerate stater, we investigate the attraction basins forr andr * , which turn out to represent the physically interesting states, namely the pure phases including the metastable states and the coexistence states. We do not consider here the hyperbolic equilibrium states, since they are unstable.
We introduce the following functions, where the index P stands for Pure phase and C for Coexistence: where f * * is the convex hull of f , as given by (22). Notice that for all r, we have G P (r) = G C (r) for ρ / ∈]ρ * 1 , ρ * 2 [, and G P (r) ≤ G C (r) for ρ ∈ [ρ * 1 , ρ * 2 ]. These two functions are candidates to be Lyapunov functions, and we have indeed the following proposition. • In the spinodal zone, that is ρ ∈]ρ − , ρ + [, G C is a Lyapunov function on the whole domain ]0, ρ[×]ρ, 3[. • In the pure stable zones, that is ρ ≤ ρ * 1 or ρ ≥ ρ * 2 , G P is a Lyapunov function on the whole domain , there are two basins of attraction: a pure phase basin Ω P where G P is a Lyapunov function, a coexistence basin Ω C = (]0, ρ[×]ρ, 3[) \ Ω P where G C is a Lyapunov function. More precisely, in the case of metastable vapour, ρ ∈ [ρ * 1 , ρ − ], we have where ρ 1 is defined by ρ 2 (ρ 1 ) =ρ 2 defined in in Proposition 2.2-(2b), and ρ 2 (ρ 1 ) is an implicit function associated to A similar characterization holds in the liquid metastable zone.
Thus we are left with positivity, which we study for the three above zones. Spinodal zone: ρ ∈]ρ − , ρ + [. By definition of f * * we have G C (r) > 0 if r = r * . Then the attraction basin corresponds to the set of states r satisfying (23), that is ρ ∈]ρ * 1 , ρ * 2 [. Notice that for any ρ in this region G P (r) < 0 in a neighbourhood ofr.
More insight on the function ρ 2 (ρ 1 ) is obtained by rewriting (36) as The partial derivatives of ϕ read Convexity properties of f ensure that µ( On the other hand, ∂ ρ2 ϕ = 0 for ρ 1 = ρ and r = (ρ 1 ,ρ 2 ) as defined above. It remains to study the sign of . Then ∂ ρ2 ϕ < 0. Thus the implicit function theorem applies and ρ 2 (ρ 1 ) = −∂ ρ1 ϕ/∂ ρ2 ϕ < 0 for ρ 1 ∈]ρ 1 , ρ[, with infinite derivative in ρ 1 and ρ. In conclusion G P is a Lyapunov function for the attraction basin defined by ]0, Figure 2. Stable equilibrium states and attraction basins for (27): pure gaseous zone (left)pure gaseous metastable zone (center) -spinodal zone (right). The blue areas refer to nonattainable states according to (23). The stable equilibrium points appear in red. The attraction basins are in white, except for the metastable basin in green.

Isothermal Euler system
In this section we propose an example of coupling between the above thermodynamical model and hydrodynamics. This toy-model allows us to explore through examples the ability of the system to catch metastable and coexistence states.
Notice that recombining the equations in ρ, ρ 1 , ρ 2 we obtain the following relaxation equation on the volume fraction α 1 : This kind of equation is classical at least for the transport operator, see for instance [3,7], the right-hand side is specific to the relaxation model we use.
In system (37), ε is a characteristic time, which compares the relaxation time to thermodynamical equilibrium with respect to the hydrodynamical time. Three regimes can be considered here. First the infinite relaxation regime, where ε = +∞, so that only the convective part of (37) is involved, the thermodynamics is too slow to affect the hydrodynamical motion. At the opposite, the instantaneous relaxation corresponds to ε = 0, so that (ρ 1 , ρ 2 ) and (α 1 , α 2 ) in (37) are given by the long time limit of the dynamical system (27). In this regime, also known as quasi-static regime, the thermodynamical relaxation is instantaneous, giving rise to a quasilinear system. In between, finite values of ε lead to an actual relaxation system.
We now address the numerical discretization. Being given a time step ∆t and a space step ∆x, we denote as usual by W n j an approximation of W (t n , x j ), where t n = n∆t and x j = j∆x, n ≥ 0, j ∈ Z. This is complemented by discrete initial data W 0 j , j ∈ Z. System (37) is discretized using a time splitting approach. First the convective part of (37) is approximated by a conservative finite volume scheme whose generic form is where F n j+1/2 is determined by a numerical flux F (W n j , W n j+1 ). We choose here the so-called HLLC numerical flux. Recall the classical HLL flux [5] for to states W L and W R reads where . The HLLC variant consists in taking into account the eigenvalue u, see e.g. [12] for a general presentation. We define the 1st and 4th components of the HLLC flux by F HLLC 1,4 (W R , W L ) = F HLL 1,4 (W R , W L ), and the 2nd and 3rd components by In a second step we treat the relaxation source terms of (37). We consider instantaneous relaxation with ε = 0. In practice one determines in which basin of attractionr n+1 j belongs to and defines r n+1 j as the corresponding equilibrium state, according to Proposition 2.4. Since the different equilibrium states belong to hyperbolicity regions of the system (37), loss of hyperbolicity is not a concern.
All the following examples are Riemann problems on [0, 1], with W 0 (x) = W L 1 x<1/2 (x) + W R 1 x>1/2 (x), W L and W R being two constant vectors. They are computed with 10000 space steps, that is ∆x = 10 −5 , the time step being computed at each time iteration by the CFL condition ∆t n = σ CF L ∆x max j∈Z ,k∈{1,2,3,4} where 0 < σ CF L < 1 is a constant which was taken equal to 0.95.

Nucleation by compression
We begin with an example of compression of a pure vapour state, for different values of the velocity u. Initial data for all the cases we consider involve constant densities on [0, 1] and opposite velocities: so that the initial perturbation remains within the pure vapour phase. We apply successively four values of u L :  surrounding the bubble, second a very thin zone of mixture (ρ 1 = ρ * 1 , ρ 2 = ρ * 2 ). It is very clearly evidenced on the the volume fraction α profile, see Figure 4 -center left, but appears as "kinks" for instance on the velocity or partial densities pictures.
If we increase the compression, with u L = 1.5, we observe that the metastable zone eventually disappears, while the tiny mixture zone is still present, see Figure 5.

Cavitation by decompression
In a somehow symmetric way, we can obtain cavitation, that is creation of a vapour bubble from a liquid initial condition. More precisely, we consider initial densities ρ = 1.9 > ρ * 2 , with ρ 1 = 1.87 and ρ 2 = 1.92. To these data we apply a velocity u L = −u R < 0. As in the compressive case, for small values values of the velocity we deal with a classical Riemann problem for the Euler system with a monotone pressure law. These results are not shown here. We focus on higher values of the velocity for which a bubble appears, namely u L = −0.4 = u R (strong decompression, Figure 6) and u L = −2 = u R (very strong decompression, Figure 7).
The situation here is more complex than in the compression case. First notice that there is no intermediate metastable state. Next, it seems that the rarefaction induced by the decompression gives rise to composite waves, see Figure 7. The tiny mixture zone exists in both cases, see Fig. 6 and 7 -center left the void fraction α 1 . The situation is more sensitive to the CFL condition: notice the overshoots in the pressure profile in Figure 7 -center left: they should not go over p * in the coexistence zone. If we reduce the CFL condition (e.g. σ CF L = 0.65), the pressure remains bounded by p * .