A phase field model for liquid-vapour phase transitions

We propose a model describing the liquid-vapour phase transition according to a phase-field approach. The model takes up a setting proposed by the second author, where a phase field is introduced whose equilibrium values 0 and 1 are associated to the liquid and vapour phases. The phase field obeys Ginzburg-Landau equation and enters the constitutive relation of the density, accounting for the sudden density jump occurring at the phase transition. In this paper we concern ourselves with the extension of this model to take into account the existence of the critical point in the coexistence line.


Introduction
In this paper an idea put forward in [1] concerning a phase field description of water-vapour phase transition is developed. The main macroscopic manifestation of the liquid vapour phase transition (and, to a smaller extent, of the ice-water transition) is a sudden jump at a given pressure in the density of the substance. The starting assumption made in the cited paper is the decomposition of the density ρ of the fluid into two contributions ρ 0 and ρ 1 , such that The contribution ρ 0 is taken as a function ρ 0 =ρ 0 (p, θ) of the pressure p and temperature θ, while ρ 1 also hinges on an independent phase field ϕ. Such a decomposition entails the attainment of the Andrews curves with their characteristic plateaus. Otherwise, these would be obtained from the van der Waals state equation supplemented by the Maxwell construction. In the cited approach, instead, it is the variation of the order parameter which gives rise to the horizontal piece of the isothermal lines (which prevents the relation between pressure and specific volume to be one-to-one). The additional phase field requires a further evolution equation (besides continuity, Navier-Stokes and heat equations) which is of course the time-dependent Ginzburg-Landau equation. We remark that in the classical Ginzburg-Landau approach to the liquid-vapour critical point, the order parameter is identified with the difference ρ−ρ c between the density ρ and the critical density ρ c (see e.g. [2,3]). Here, the order parameter is independent of the density of the fluid and stands out from the constitutive decomposition (1) of the overall fluid density.
In this paper we concern ourselves with the modeling of the liquid-vapour phase transition along the same lines of [1], but also including the peculiarities due to the presence of the critical point. This aspect was not fully considered in the previous model, which makes sense specifically for temperatures below the critical temperature. As a consequence of the existence of the critical point, there always exists a sequence of equilibrium states which connect in a continuous way the liquid and the vapour phase (see Fig. 1). Taking into account this kind of transformations, it is not possible to assign a unique fixed value of ϕ to the liquid and vapour phases anymore. We remark that, due to the existence of the critical point, it is possible to speak of distinct phases only on the coexistence curve. So it is necessary to consider a phase field which admits a continuous set of equilibrium values under different temperature and pressure conditions. Obviously this exhibits the second order (or continuous) character of the transition at the critical point. In usual phase field models (see e.g. [1,4]) the homogeneous phases of a substance are associated to fixed discrete values of the order parameter (for example ϕ = 0 and ϕ = 1). Intermediate values are associated to the thin interfacial layers between two phases; in particular, they are associated, in stationary conditions, with non vanishing values of ∇ϕ. This will not be true anymore if we assume a phase field which have a continuous set of equilibrium values associated to homogeneous phases, similarly to a Ginzburg-Landau order parameter in a strict sense (for second order phase transition). In this case intermediate values of the phase field can be associated both to interfacial layers and to homogeneous phases, close to the critical point.

Model equations and thermodynamical consistence
We start by setting the general thermodynamic framework of the model we are going to propose. Our description is based on the introduction of an order parameter ϕ connected with the density jump of the transition. The equation for ϕ is formulated as an additional balance equation to which is associated a corresponding balance of powers [1,5,6]. As a consequence, a variation of ρ 1 is associated to a supplementary power flux in addition to the classical mechanical power pρ/ρ.
We construct an evolutive model based on four balance equations for the fundamental thermomechanical fields of the fluid, that is: the phase field ϕ, the density ρ, the macroscopic velocity v and the temperature θ. The thermomechanical state of the system is defined by the set of variables where h = ∇ϕ. We remark that we are considering a compressible fluid, both in its liquid and in its vapour phase. In this situation, the choices of the pressure or the density ρ as state variables are equivalent. From our point of view, it is more natural to choose the pressure as state variable, because it appears explicitly in the Ginzburg-Landau equation for the order parameter. We state the four balance equations of the model in the order: order balance equation, mass balance, momentum balance and energy balance where: -the superimposed dot represents the material derivative of a field:φ = ∂φ/∂t + v · ∇φ; -k is the internal order production and p the order flux. For the meaning of this terminology we refer to [5]. The relation of these fields with the order parameter ϕ is given in Eq. (9); -T E is the extra stress (see Eq. (14)), and b the external force density; -e is the internal energy density, P i is the internal power density which is the sum of the internal mechanical power P i m associated to motion equation and the power P i ϕ associated to order balance; h is the internal thermal power density, which obeys its own balance in terms of heat flux q and heat supply r: The balance laws have to be completed by appropriate constitutive relations for k, p, T E , ρ. The phase field equation suitable for the description of a phase transition takes the form of a Ginzburg-Landau equation: where τ and κ are positive parameters, f is the Ginzburg-Landau potential and f ϕ = ∂f /∂ϕ. Hence the fields k and p satisfy the constitutive equations Equation (8) allows us to write the balance of powers associated to order parameter by multiplying both members forφ and extracting a divergence contribution: We identify the first member with the internal power of the order parameter; using the identity ∇φ =ḣ + h · ∇v.
we find useful to write the internal power as where D = 1 2 (∇v + ∇v T ) is the rate of deformation tensor. The internal mechanical power P i m has the classical expression and the extra-stress is assumed in the form It includes a viscous contribution and a coupling term with the order parameter gradient h which is necessary to assure the thermodynamical consistence (see below). The viscosity µ(θ, ϕ) is a function of the temperature and the phase which vanishes in the vapour and gas phases. The constitutive equation for ρ is conveniently expressed in terms of the specific volume ν := 1/ρ. We will assume a generic constitutive relation ν = ν(p, θ, ϕ) 1 . The choice of the functionν is strictly related to the function f (p, θ, ϕ) by the thermodynamical constraints and will be specified after we will have examined the consequences of the Second Law. We find it convenient to express the mechanical power (13) in terms of the specific volume ν Now we examine the constraints on the constitutive equations due to the Second Law, in the form of the Clausius-Duhem inequality. In terms of the free energy ψ = e − θη, the inequality can be expressed as where g = ∇θ. Substituting the expressions (12) and (15) we obtain Firstly we observe that constitutive equation (14) makes the term proportional to D a non-negative definite term, provided that µ ≥ 0. From inequality (17), it is plain the reason to include the contribution proportional to h ⊗ h in the extra stress tensor. In order to satisfy the Clausius-Duhem inequality in the most simple (non trivial) way, we add the following constitutive relations Then, the following expressions are easily obtained 2 takes the meaning of a specific Gibbs' free energy and ψ = Φ − pΦ p . Substituting the expressions ofν and T E in the internal power expression we obtain Then, the first law of thermodynamics gives us the heat equation We resume the set of equations describing the model The last equation is the constitutive equation for the density. It could be regarded as constitutive equation for the pressure if we would consider the density as a state variable. The Clausius-Duhem inequality turns out to be very useful in establishing the coupling between the different equations. It has nothing to say about the choice of the function f instead.

Equilibrium properties.
From the Ginzburg-Landau equation (8), the equilibrium valueφ =φ(p, θ) of the order parameter is a solution of the equation which also should be a local minimum of f , for stability reasons. As the system has at most two coexisting phases, f (p, θ, ϕ) has to admit, as a function of ϕ, at most two local minimaφ 1 ,φ 2 . On the coexistence line the two local minima have to be absolute minima, that is: In general, there is a strip-like region bordering the coexistence line in which both minima are present: one of them is an absolute minimum and represents the true equilibrium state in that region, while the other one is a relative minimum and represents a metastable state, which can be a superheated fluid or an undercooled vapour. Elsewhere, the order parameter has a unique equilibrium value, in particular for θ > θ c (or p > p c ). The equilibrium Gibbs free energy (26) for the homogeneous substance (that is for h = 0) in the i-phase isΦ and the coexistence line is characterized byΦ 1 =Φ 2 . The usual thermodynamics of phase equilibrium is then recovered. Firstly, we remark that, thanks to (31), we have So, the equilibrium equation of state of the i-phase which follows from (24) is represented by the usual relation Similarly, from (25), the equilibrium entropy is given bȳ The specific volume jump due to phase change is one of the most important parameters of the transition, together with the latent heat These quantities are related by the well known Clausius-Clapeyron equation which can be easily obtained deriving both members of equilibrium condition with respect to θ and using (36), (37).

The Ginzburg-Landau potential
We now turn to the choice of the Ginzburg-Landau potential f . The phase field approach to first order phase transitions involves an order parameter ϕ which assumes a fixed value (say ϕ = 0) in one of the two phases, and a different fixed value (say ϕ = 1) in the other phase. The order parameter does not suffer jump discontinuities at phase boundaries, but varies smoothly across a interfacial layer from one equilibrium value to the other [6,7,8]. This is obtained by subjecting the order parameter to a Ginzburg-Landau equation, with a potential f (ϕ) having minima only at ϕ = 0, 1 for all values of the temperature and pressure (or any other relevant physical field). In the paper [1] it was chosen a potential such that it had always two relative minima at ϕ = 0 and ϕ = 1. That minima were associated respectively to the liquid and vapour phases. This setting is not suited to model the phases in the neighborhood of the critical temperature. In fact, at the critical temperature θ c , the liquid and vapour phases merge in a unique phase which is sometimes named gas phase in a strict sense 3 (or supercritical fluid in the region θ > θ c , p > p c ). More precisely, at the critical point, a second order phase transition occurs. We cannot associate a fixed value of the order parameter to the states lying on one side of the coexistence line, if this one have an end point. In fact, consider the transformation represented in Fig. 1. No phase transition occurs, so the order parameter has no jump discontinuity along that way. Therefore the order parameter should have a continuous set of equilibrium values. In particular, this rules out an approach based on three phases, liquid, vapour and gas, associated to fixed values of the order parameter (for example, 1, -1, 0): in that case, a physically inexistent first order phase transition between gas and liquid and between gas and vapour would take place somewhere in the phase diagram. The most direct generalization of the approach followed in [1] is to assume a Ginzburg-Landau potential with minimal points varying continuously with temperature at θ < θ c . The two minima which are present in the coexistence region of liquid and vapour merge in a unique minimum ϕ = 0 for θ ≥ θ c . The pressure will not move the position of the minima, but can make the order parameter jump from the liquid minimum to the vapour minimum and vice versa. So, in isothermal condition and inside the one phase region, the order parameter will be effectively constant and the Navier-Stokes equation is in fact decoupled from Ginzburg-Landau equation. We start by considering the double well potential having two minima at x = ±u and a maximum at x = 0. We next define the C 1 piecewise polynomial and odd function Any linear combination F (x; u) + hG(x; u), with a real coefficient h, is a function of x having a minimum or a flexus in x = ±u and no other minimum.
The Ginzburg-Landau free energy density is taken in the form whereĥ is a function vanishing on the coexistence line, namelyĥ(p 0 (θ), θ) = 0. Then, the specific volume is given by We identify the components of the decomposition (1) in the following way Then, the difference in volume on the coexistence line, corresponding to the jumping in ϕ from the minimumφ 1 = u to the minimumφ 2 = −u, is If we assumeĥ(p, θ) = p − p 0 (θ) for θ < θ c , then the parameter u is directly related to the measurable volume differencê The liquid phase corresponds to the phase diagram region θ < θ c , p > p 0 (θ), that is, in the case we are considering, to h > 0. For h > 0 the absolute minimum is attained at ϕ = u, so, for this choice ofĥ, the liquid phase is associated to the right minimum, while the vapour phase is associated to the left minimum. For θ > θ c we have ∆ν = 0 and soû(θ) = 0. We note that G(x; 0) = 0, so, for θ > θ c the Ginzburg-Landau potential is which has the only minimum ϕ = 0. We associate that value to the gas phase. On isothermal and homogeneous conditions (so that ∇ϕ = 0) and with externally prescribed homogeneous pressure field p(t), the evolution of the volume is ruled by the equations For τ 1, one can assume that at every instant the order parameter is equal to the equilibrium valueφ(t) = ±u, so that the equilibrium equation of state is ν(p, θ) = (f 0 ) p (p, θ) +ĥ p (p, θ)G(φ; u).
(52) u h h=u Figure 2: Structure of the minima of the potential f (p, θ, ϕ) at varying h > 0, u > 0. In the region h < 0 the structure is the same with the interchange of ϕ with −ϕ.
The equilibrium valueφ considered as a function of p, θ, is not a single-valued function, since there is an hysteresis behaviour in the interval |h| < u. Considering only the stable states and discarding the metastable ones, amounts to takeφ asφ Assuming this functionφ(p, θ) the equilibrium isotherm at θ < θ c in the ν − p plane is which has a flat piece at pressure p 0 (θ) with length ∆ν(θ) = 4u 3 /3. For pressures p = p 0 (θ) the curve is determined essentially by f 0 (p, θ). For θ ≥ θ c ,φ = 0 and ν = (f 0 ) p .
A difficulty in this model is given by a possible diverging contribution in the equilibrium entropy if ∆ν(θ) is singular at θ = θ c , as, for example, if ∆ν ∼ |θ − θ c | β , with β < 1. In fact, consider the expression of the entropy η(p, θ, ϕ) = −f θ (p, θ, ϕ) = η 0 (p, θ) + uu θ ϕ 2 − h θ G(ϕ; u) + 2huu θ ϕ, where η 0 (p, θ) := ∂f0 ∂θ . So, a diverging expression for u θ would give rise to a diverging entropy at critical point. We can directly relate the equilibrium entropyη with ∆ν(θ). At the equilibrium values ϕ = ±u, we have the entropies or, remembering that ∆ν(θ) = 4u 3 /3, In particular, at the critical point u = 0, ∆ν = 0 and the entropy is given bȳ Now, in general (∆ν) θ is a singular quantity at the critical point; for example, according to the Landau theory of the second order phase transitions, ∆ν ∝ |θ − θ c | 1/2 . So, in this case we would obtain a divergent entropy at the critical temperature. In paper [1] both equilibrium and non-equilibrium entropy function was a regular function of the temperature; nevertheless that model involved a temperature dependence of the volume jump near critical point in the form ∆ν ∼ |θ − θ c |. Furthermore, the latent heat does not approach to zero at critical point. This shortcoming can be avoided in ad hoc way assuming a cancelling divergence contribution in η 0 or following more closely the Landau theory of second order phase transitions, in which the entropy remains a regular function at the critical point, while the equilibrium value of the order parameter has a square root singularity at critical point. We suggest to consider the potential with a > 0 and u(θ) an increasing function of θ such that as, for example, in the following class of functions: , q > 0, β > 0.
As before, h has to satisfy h(p 0 (θ), θ) = 0 and h ≶ 0 for p ≶ p 0 (θ). The minima are solutions of the third order algebraic equation The same algebraic equation for the minima would be obtained (except for solutionsφ ≥ ±1) from the fourth order polynomial potential which is the classical form of the Ginzburg-Landau potential for second order phase transitions. The employment of the logarithmic potential (also used in different contexts, see for example [9]) has the advantage that the only attainable minima are constrained in the interval (−1, 1). This would not be of great importance for the only study of critical region, where order parameter is small: ϕ 1. On the contrary, if one wants the phase field to be bounded under all the temperature and pressure conditions, the fourth order polynomial potential would not be suitable. If u < 0, there are two local minima −1 <φ Figure 3: Structure of the minima of the potential f (p, θ, ϕ) at varying h > 0, u > 0. In the region h < 0 the structure is the same with the interchange of ϕ with −ϕ. ϕ 2 < 1 for |h/a| <h(u), while there is only one local minimum −1 <φ < 1 if |h/a| ≥h(u). (The functionh(u) has not a simple expression, but 2|u/3| 3/2 < h < |u|; see Fig. 3). For u > 0, there is always only one minimum in [−1, 1]. In particular, consider the case h = 0; for θ < θ c (on the coexistence line), the minima are given byφ while for θ > θ c (on the critical isochore) the only minimum isφ = 0. For h = 0, the minima have complicated algebraic expressions which are not very illuminating; we only remark that the absolute minimum satisfies the following estimate: In Fig. 4 it is plotted the value ofφ in function of h for three values of u. The which, owing to algebraic equation (62) for the equilibrium valueφ, implies the following state equation The specific volume jump at the transition is ∆ν = 4h p |u| 1/2 .
If u ∼ (θ/θ c −1), the quantity ∆ν has the typical mean-field singularity |θ−θ c | 1/2 near the critical point, needlessly to have u(θ) also singular. Finally, equilibrium entropy is given bȳ Now, u(θ) is a regular function of θ, and soη is not diverging in the neighbourhood of the critical point.
With this choice of the Ginzburg-Landau potential, the equilibrium value (absolute minimum) of the order parameterφ is not constant in the inner of the regions of the phase diagram marked as 'liquid' and 'vapour'. In this sense, the Ginzburg-Landau equation is never completely decoupled from the Navier-Stokes equation. However, as the order parameter is approximatively ±1 in the region |h| a (far from transition line) and/or |u| ∼ 1 (well below the critical temperature).