Singularities in Euler flows: multivalued solutions, shock waves, and phase transitions

In this paper, we analyze various types of critical phenomena in one-dimensional gas flows described by Euler equations. We give a geometrical interpretation of thermodynamics with a special emphasis on phase transitions. We use ideas from the geometrical theory of PDEs, in particular, symmetries and differential constraints to find solutions to the Euler system. Solutions obtained are multivalued, have singularities of projection to the plane of independent variables. We analyze the propagation of the shock wave front along with phase transitions.


Introduction
Various types of critical phenomena, such as singularities, discontinuities, wave fronts and phase transitions, have always been of interest from both mathematical [1,2,3] and practical [4] viewpoints. In the context of gases, discontinuous solutions to the Euler system, describing their motion, are usually treated as shock waves. In the past decades such phenomena have widely been studied, see for instance, [5] for the case of Chaplygin gases, [6,7], where the weak shocks are considered. It is worth mentioning also [8,9], where the influence of turbulence on shocks and detonations is emphasized.
This paper can be seen as a natural continuation of [10], where we have considered the case of ideal gas flows. Here, we elaborate the case of more complicated and at the same time more interesting from the singularity theory viewpoint model, the van der Waals model, known as one of the most popular in the description of phase transitions. So, singularities of shock wave type that can be viewed as in some sense singular solutions to the Euler system are analyzed together with singularities of purely thermodynamic nature, phase transitions. Our approach to finding and investigating such phenomena is essentially based on the geometric theory of PDEs [11,13,12,15,14]. Namely, we find a class of multivalued solutions to the Euler system (see also [16]), and singularities of their projection to the plane of independent variables are exactly what drives the appearance of the shock wave [17]. Similar ideas were used in a series of works [18,19,20], where multivalued solutions to filtration equations were obtained along with analysis of shocks. To find such solutions we use the idea of adding a differential constraint to the original PDE in such a way that the resulting overdetermined system of PDEs is compatible [21]. The same concepts were also used in [22], where a general solution to the Hunter-Saxton equation was found, in [23], where the two-dimensional Euler system was considered, and in [24], where this approach was applied to the Khokhlov-Zabolotskaya equation.
The paper is organized as follows. Section 2 is preliminary, we describe there the necessary concepts from thermodynamics. In Section 3, we analyze a multivalued solution to Euler equations and its singularities, including shock waves and phase transitions. In the last section, we discuss the results. The essential computations for this paper were made with the DifferentialGeometry package [25] in Maple.

Thermodynamics
In this section, we give necessary concepts from thermodynamics. As we shall see, geometrical interpretation of thermodynamic states allows one to use Arnold's ideas from the theory of Legendrian and Lagrangian singularities [1,2,3], which are crucial in description of phase transitions. The geometrical approach to thermodynamics was initiated already by Gibbs [26], and was further developed in, for example, [27,28], and recently in [29]. For more detailed analysis we also refer to [30].

Legendrian and Lagrangian manifolds
Consider the contact space R 5 , θ with coordinates (s, e, ρ, p, T ) standing for specific entropy, specific inner energy, density, pressure and temperature. The contact structure θ is given by Then, a thermodynamic state is a Legendrian manifold L ⊂ R 5 , θ , i.e. θ L = 0 and dim L = 2. From the physical viewpoint this means that the first law of thermodynamics holds on L. Due to (1), it is natural to choose (e, ρ) as coordinates on L. Then, a two-dimensional manifold L ⊂ R 5 , θ is given by where the function S(e, ρ) specifies the dependence of the specific entropy on e and ρ. Note that determining a Legendrian manifold L by means of (2) requires the knowledge of S(e, ρ), while in experiments one usually gets relations between pressure, density and temperature. To this reason, we get rid of the specific entropy s by means of projection π : R 5 → R 4 , π(s, e, ρ, p, T ) = (e, ρ, p, T ) and consider an immersed Lagrangian manifold π L = L ⊂ R 4 , Ω in a symplectic space R 4 , Ω , where the structure symplectic form Ω is Then, one can treat thermodynamic state manifolds as Lagrangian manifolds L ⊂ R 4 , Ω , i.e. Ω |L = 0. In coordinates (T, ρ), a thermodynamic Lagrangian manifold L is given by two functions Since Ω |L = 0, the functions P (T, ρ) and E(T, ρ) are not arbitrary, but are related by where [f, g] is a Poisson bracket of functions f and g on R 4 , Ω uniquely defined by the relation Equation (4) forces the following relation between P (T, ρ) and E(T, ρ): (−ρ −2 T −1 P )T = (T −2 E)ρ, and therefore the following theorem is valid: The Lagrangian manifold L is given by means of the Massieu-Planck potential φ(ρ, T ) Remark 2 Having given the Lagrangian manifold L by means of (3), one can find out the entropy function S(e, ρ) solving the overdetermined system Se with compatibility condition (4).

Riemannian structures, singularities, phase transitions
There is one more important structure arising, as it was shown in [29], from measurement approach to thermodynamics. Indeed, if one considers equilibrium thermodynamics as a theory of measurement of random vectors, whose components are inner energy and volume v = ρ −1 , one drives to the universal quadratic form on (R 4 , Ω) of signature (2, 2): where · is the symmetric product, and areas on L, where the restriction κ|L of κ to L is negative, are those where the variance of a random vector (e, v = ρ −1 ) is positive [29,31]. Using (5), we get and taking into account (5), we conclude that the condition of positive variance is satisfied at points on L, where eT > 0, pρ > 0, which is known as the condition of the thermodynamic stability. Let us now explore singularities of Lagrangian manifolds. We will be interested in the singularities of their projection to the plane of intensive variables (p, T ), i.e. points where the form dp ∧ dT degenerates. We will assume that extensive variables (e, ρ) may serve as global coordinates on L, i.e. the form de ∧ dρ is non-degenerated everywhere. The set where dp ∧ dT = 0 coincides with that where 2ρ −1 φρ + φρρ = 0, or, equivalently, where the from κ|L degenerates. A manifold L turns out to be divided into submanifolds Li, where both (e, ρ) and (p, T ) may serve as coordinates, or, equivalently, the form (6) is non-degenerated. Such Li are called phases. Additionally, those of Li, where (6) is negative, are called applicable phases. Thus we end up with the observation that singularities of projection of thermodynamic Lagrangian manifolds are related with the theory of phase transitions. Indeed, by a phase transition of the first order we mean a jump from one applicable state to another, governed by the conservation of intensive variables p, T , and specific Gibbs potential γ = e − T s + p/ρ, which in terms of the Massieu-Planck potential is expressed as γ = −T (φ+ρφρ) [30]. Consequently, to find the points of phase transition, one needs to solve the system where p and T are the pressure and the temperature of the phase transition, and ρ1 and ρ2 are the densities of gas and liquid phases.

Example 3 (Ideal gas)
The simplest example of a gas is an ideal gas model. In this case the Legendrian manifold is given by where R is the universal gas constant, and n is the degree of freedom. The differential quadratic form κ|L is It is negative definite on the entire L, and there are no phase transitions, nor are there singularities of projection of L to the p − T plane.
Example 4 (van der Waals gas) To define the Legendrian manifold for van der Waals gases we will use reduced state equations: The differential quadratic form κ|L is In this case it changes its sign, the manifold L has a singularity of cusp type. The singular set of L, called also caustic, and the curve of phase transition are shown in Figure 1.

Euler equations
In this paper we will study non-stationary, one-dimensional flows of gases, described by the following system of differential equations: • Conservation of momentum • Conservation of mass • Conservation of entropy along the flow Here u(t, x) is the flow velocity, ρ(t, x) is the density of the medium, and s(t, x) is the specific entropy. System (10)- (12) is incomplete. It becomes complete once extended by equations of thermodynamic state (2). We will be interested in homentropic flows, i.e. those with s(t, x) = s0. On the one hand, this assumption satisfies (12) identically, on the other hand, it allows us to express all the thermodynamic variables in terms of ρ. Indeed, the entropy s has the following expression in terms of the Massieu-Planck potential φ(T, ρ): s = φ + T φT [30]. Putting s = s0, we get an equation s0 = φ + T φT , which determines T (ρ) uniquely, since its derivative with respect to T is positive due to the negativity of κ|L. Substituting T (ρ) into (3), one gets p = p(ρ). Thus we end up with the following two-component system of PDEs: where A(ρ) = p ′ (ρ)/ρ. We do not specify the function A(ρ) yet, we shall do this while solving (13).

Finding solutions
To find solutions to system (13), we use the idea of adding a differential constraint to (13), compatible with the original system. It is worth mentioning that a solution is an integral manifold of the Cartan distribution on (13) (see [11,13,12] for details). This geometrical interpretation of a solution to a PDE allows to find ones in the form of manifolds, which, in general, may not be globally given by functions. This approach gives rise to investigation of singularities in a purely geometrical manner, which is shown in this paper. In general, finding differential constraints is not a trivial problem. But having found ones, the problem of finding solutions is reduced to the integration of a completely integrable Cartan distribution of the resulting compatible overdetermined system. In case the Cartan distribution has a solvable transversal symmetry algebra, whose dimension equals the codimension of the Cartan distribution, we are able to get explicit solutions in quadratures by applying the Lie-Bianchi theorem (for details we refer to [11,13,12]).
We will look for a differential constraint compatible with (13) in the form of a quasilinear equation where functions α(ρ) and β(ρ) are to be determined. We will denote system (13)-(14) by E .
The proof of Theorem 5 is more technical rather than conceptual. First of all, we lift system (13)- (14) to the space of 3-jets J 3 (R 2 ) by applying total derivatives Dt = ∂t + ut∂u + ρt∂ρ + utt∂u t + ρtt∂ρ t + . . . , to equations of E the required number of times consequently. The resulting system E3 ⊂ J 3 (R 2 ), consisting of equations only of the third order, contains 9 equations for 8 variables of purely third order, i.e. uttt, uxxx, utxx and so on. Eliminating them from E3, we get 7 relations (6 obtained by lifting E to J 2 (R 2 ) plus one remained from eliminations of third-order variables). Again, we eliminate all the variables of the second order and we get four relations of the first order. Eliminating ux, ut and ρt we end up with the expression of the form ρ 3 x G(ρ, u) = 0, where G(ρ, u) is a polynomial in u, whose coefficients are ODEs on α(ρ), β(ρ) and A(ρ), solving which we get (15). It is worth saying that these computations are algebraic and are well suited for computer algebra systems. (8) and (9), one can show that the function A(ρ) = p ′ (ρ)/ρ given in (15) corresponds to that of • ideal gas in case of

Remark 6 Using
• van der Waals gas in case of The case of ideal gases was thoroughly investigated in [10]. Here, we are interested in the case of van der Waals gases.
Summarizing, we have a compatible overdetermined system of PDEs where functions α(ρ), β(ρ) and A(ρ) are specified in (15). This system is a smooth manifold E in the space of 1-jets J 1 (R 2 ) of functions on R 2 . Since dim J 1 (R 2 ) = 8, and E consists of 3 relations on J 1 (R 2 ), dim E = 5. The dimension of the Cartan distribution CE on E equals 2, therefore codim CE = 3. Let us choose (t, x, u, ρ, ρx) as internal coordinates on E . Then the Cartan distribution CE is generated by differential 1-forms where ρxx, ρxt, ut, ux, ρt are expressed due to E and its prolongation where A(ρ) is given by (15). We shall look for integrals of the distribution (17)- (22), which give us an (implicit) solution to (13)- (14).
Let us choose another basis κ1, κ2, κ3 in CE by the following way: Due to the structure of the symmetry Lie algebra g, the form κ1 is closed [11,12], and therefore locally exact, i.e. κ1 = dQ1, where Q1 ∈ C ∞ (J 1 ), while restrictions κ2|M 1 and κ3|M 1 to the manifold M1 = {Q1 = const} are closed and locally exact too. Integrating the differential 1-form κ1 we observe that variables u, ρ, t, x can be chosen as local coordinates on M1 and where α1 is a constant. Integrating restrictions κ2|M 1 and κ3|M 1 , we get two more relations that give us a solution to (13)- (14) implicitly: and where we have already substituted A(ρ) from (15), and α2, α3 are constants. The graph of a multivalued solution for the density is shown in Figure 2. We used substitution (16), where C5 = 240, n = 3, together with C2 = 1, α1 = 1, α2 = 2, α3 = 1.

Phase transitions
Having a solution, one can remove the phase transition curve from the space of thermodynamic variables to R 2 (t, x). Indeed, on the one hand, we have all the thermodynamic parameters as functions of (t, x). On the other hand, we have conditions on phase transitions (7) in the space of thermodynamic variables. In combination they give us a curve of phase transitions in (t, x) plane. Phase transitions together with the shock wave are presented in Figure 4.

Discussion
In the present work we analyzed critical phenomena in gas flows of purely thermodynamic nature, which are phase transitions, and shock waves arising from singularities of solutions to the Euler system. To obtain such solutions we used a differential constraint compatible with the original system. In this work it was found in a purely computational way, and it seems interesting how to get it in a more regular way. One of possible ways to find such constraints is using differential invariants. Then, constraints can be found constructively by solving quotient PDEs, which was successfully realized in [22]. We hope to make use of this method in the future research. The analysis of phase transitions showed that sometimes shock waves can be accompanied with phase transitions, which is shown in Figure 4, since the phase transition curve intersects the shock wave front, and on the one side of the discontinuity curve we observe a pure gas phase, while on the other side we can see a wet steam.