Well-posedness of an extended model for water-ice phase transitions

We propose an improved model explaining the occurrence of high stresses due to the difference in specific volumes during phase transitions between water and ice. The unknowns of the resulting evolution problem are the absolute temperature, the volume increment, and the liquid fraction. The main novelty here consists in including the dependence of the specific heat and of the speed of sound upon the phase. These additional nonlinearities bring new mathematical difficulties which require new estimation techniques based on Moser iteration. We establish the existence of a global solution to the corresponding initial-boundary value problem, as well as lower and upper bounds for the absolute temperature. Assuming constant heat conductivity, we also prove uniqueness and continuous data dependence of the solution.


Introduction
In the present contribution we prove the well-posedness of an initial-boundary value problem associated with the following system coupling a quasi-linear parabolic internal energy balance (for the absolute temperature θ ) with an integro-differential equation for the relative volume increment U , and a differential inclusion ruling the evolution of the phase variable χ as follows: c(χ)e 1 (θ) t − div (κ(χ)∇θ) = c ′ (χ)χ t (f 1 (θ) − e 1 (θ)) with U Ω (t) = Ω U (x, t) dx. In the previous paper [12] we have already given a motivation and a complete study of equilibria for this system, which models the water freezing in an elastic container, taking into account differences in the specific volume, specific heat and speed of sound in the solid and liquid phases. The derivation of the system from physical principles and the meaning of the symbols will be explained below in the next Section 2.
Here, we describe the mathematical difficulties and comment on previous results related to this type of systems.
There is an abundant classical literature on phase transition processes, see e.g. the monographs [2], [5], [14] and the references therein. It seems, however, that only few publications take into account different mass densities/specific volumes of the phases. In [6], the authors proposed to interpret a phase transition process in terms of a balance equation for macroscopic motions, and to include the possibility of voids. Well-posedness of an initial-boundary value problem associated with the resulting PDE system is proved there and the case of two different densities ̺ 1 and ̺ 2 for the two substances undergoing phase transitions has been pursued in [7].
Here, we deal exclusively with physically measurable quantities. All parameters have a clear physical meaning and the derivation is carried out under the assumption that the displacements are small. This enables us to state the system in Lagrangian coordinates (cf. [7] for a different approach to the subject).
The present model has been previously studied in [10] and [11] under the assumption that the speed of sound and the specific heat are the same in solid and in liquid. In terms of the system (1.1)- (1.3), this corresponds to choose constant functions λ(χ) ≡ λ and c(χ) ≡ c. For this particular case, we have proved in [10] and [11] the existence and uniqueness of global solutions, as well as the convergence of the solutions to equilibria. In reality, the specific heat in water is about the double, while the speed of sound in water is less than one half of the one in ice. The main goal of this contribution is to give a well-posedness result for a boundary value problem associated with (1.1)-(1.3) including these dependences into the model. The main result is stated in Section 4. The dependence of speed of sound and of the specific heat on the phase is expressed in terms of additional nonlinearities in the equations which have to be suitably handled. Moreover, here we also generalize the results of [10] and [11] allowing for non constant external pressure and temperature. Finally, we proceed here with a different technique for the proof of existence of solutions with respect to [10] and [11]. Since the contraction argument does not work in our situation, we discretize in time our problem (cf. Subsection 5.1), preparing thus necessary tools for future numerical investigations on this model, and prove the convergence of the scheme. The uniqueness and continuous dependence of solution on the data is proved in Section 6 following the idea already exploited in [4] where we deal with a quasi-linear internal energy balance equation coupled with a vectorial and nonlocal phase dynamic. The main estimates are obtained here by means of the energy inequality which still holds true at the discrete level (cf. Subsection 5.3). Finally, it is worth noting that a time dependent positive lower bound for the θ -component of the solution independent of the time step is established on the time discrete approximation in Subsection 5.2, while we obtain a uniform in time upper bound on the solution θ by means of a proper Moser estimate (cf. Subsection 5.4).

Balance equations
Referring to [12] for the complete deduction of the model, we consider a liquid substance contained in a bounded connected container Ω ⊂ R 3 with boundary of class C 1,1 . The state variables are the absolute temperature θ > 0, the displacement u ∈ R 3 , and the phase variable χ ∈ [0, 1]. The value χ = 0 means solid, χ = 1 means liquid, χ ∈ (0, 1) is a mixture of the two.
We make the following modeling hypotheses.
(A1) The displacements are small. Therefore, we state the problem in Lagrangian coordinates, in which mass conservation is equivalent to the condition of a constant mass density ̺ 0 > 0.
(A2) The substance is isotropic and compressible; the speed of sound and the specific heat may depend on the phase χ .
(A3) The evolution is slow, and we neglect shear viscosity and inertia effects.
(A5) The liquid phase is the reference state, and the specific volume V i of the solid phase is larger than the specific volume V w of the liquid phase.
We thus consider the evolution system consisting of a the mechanical equilibrium equation (2.1), energy conservation law (2.2), and a phase dynamic equation (2.3), where the coefficient γ 0 determines the speed of the phase transition. By (A4), the stress has the form σ = −p δ and the scalar quantity is the pressure. Here ν > 0 is a volume viscosity coefficient, λ(χ) is the Lamé constant, which may depend on χ by virtue of (A2), α = (V i −V w )/V w is a positive phase expansion coefficient by (A5), while β is the thermal expansion coefficient, which is assumed to be constant, and f vol is a given volume force density (the gravity force) with standard gravity g and vector δ 3 = (0, 0, 1). We denote by e the specific internal energy, s is the specific entropy, and q is the heat flux vector that we assume for simplicity in the form with heat conductivity κ(χ) > 0 depending possibly on χ . We assume the specific heat c V (χ, θ) in the form This is still a rough simplification, and further generalizations are desirable. According to [9,Chapter VI] or [13,Section 5], the purely caloric parts e cal and s cal of the specific internal energy and specific entropy are given by the formulas e cal (χ, θ) = c 0 (χ)e 1 (θ), s cal (χ, θ) = c 0 (χ)s 1 (θ), with Then, the specific free energy f = e − θs satisfies the conditions σ e = ̺ 0 ∂ ε f , s = −∂ θ f . With a prescribed constant latent heat L 0 and freezing point θ c > 0 at standard atmospheric pressure P stand , the specific free energy f necessarily has the form andf is a arbitrary function of χ (integration "constant" with respect to θ and ε ). We choosef so as to ensure that the values of χ remain in the interval [0, 1], and that the phase transition under standard pressure takes place at temperature θ c . More specifically, we setf where I is the indicator function of the interval [0, 1]. For specific entropy s and specific internal energy e we obtain The equation for the phase χ is obtained by assuming that −χ t is proportional to ∂ χ f with proportionality coefficient (relaxation time) γ 0 (θ) > 0, where ∂ χ is the partial Clarke subdifferential with respect to χ .
Then, the equilibrium equation (2.1) can be rewritten in the form ∇p = f vol , hence, as Ω is connected, where P is a function of time only, which is to be determined. Recall that in the reference state ε : δ = ε t : δ = 0, χ = 1, and at standard pressure P stand , the freezing temperature is θ c . We thus see from (2.4) that P (t) is in fact the deviation from the standard pressure. We assume also the external pressure in the form P ext = P stand + p 0 with a given deviation p 0 (x, t). The normal force acting on the boundary is (P (t)−̺ 0 g x 3 −p 0 )n, where n denotes the unit outward normal vector. We assume an elastic response of the boundary, and a heat transfer proportional to the inner and outer temperature difference.
On ∂Ω , we thus prescribe boundary conditions for u and θ in the form with a given symmetric positive definite matrix k(x) (elasticity of the boundary), positive functions h(x) (heat transfer coefficient), and θ Γ (x, t) > 0 (external temperature). This enables us to find an explicit relation between div u and P . Indeed, on ∂Ω we have by (2.14) that u · n = (P (t) − ̺ 0 g x 3 − p 0 (x, t))k −1 (x)n(x) · n(x). Assuming that k −1 n · n belongs to L 1 (∂Ω), we set and obtain by Gauss' Theorem that where P 0 (t) = ∂Ω p 0 (x, t)k −1 (x)n(x) · n(x) dσ(x). Under the small strain hypothesis, the function div u describes the local relative volume increment. Hence, Eq. (2.17) establishes a linear relation between the total relative volume increment U Ω (t) and the relative pressure P (t) − p 0 (x, t). We have ε : δ = div u, and thus the mechanical equilibrium equation (2.13), due to (2.4) and (2.17), reads As a consequence of (2.6), (2.9), and (2.12), the energy balance and the phase relaxation equation in (2.2)-(2.3) have the form For simplicity, we now set Note that mathematically, the subdifferential ∂I(χ) is the same as ̺ 0 L 0 ∂I(χ). The system thus reduces to the system (1.1)-(1.3) of three scalar equations -one PDE and two "ODEs" for three unknown functions θ, χ , and U , with boundary condition (2.15), (2.6).
Assuming that a solution to (1.1)-(1.3) is known with U ∈ L 2 (Ω × (0, T )), we find the vector function u by defining first Φ to be the solution to the Poisson equation ∆Φ = U with the Neumann boundary condition ∇Φ·n With this Φ , we findũ as a solution to the problem and set u =ũ + ∇Φ . Then u satisfies a.e. in Ω the equation div u = U , together with the boundary condition (2.14), that is, For the solution to (2.22)-(2.23), we refer to [8, Lemma 2.2] which states that for each g ∈ H 1/2 (∂Ω) 3 satisfying ∂Ω g · n dσ(x) = 0 there exists a functionũ ∈ H 1 (Ω) 3 , unique up to an additive function v from the set V of divergence-free H 1 (Ω) functions vanishing on ∂Ω , such that divũ = 0 in Ω ,ũ = g on ∂Ω . In terms of the system (2.22)-(2.23), it suffices to set g holds with some constants C,C . The required regularity is available here by virtue of the assumption that Ω is of class C 1,1 , provided k −1 belongs to H 1/2 (∂Ω). Note that a weaker formulation of problem (2.22)-(2.23) can be found in [1,Section 4]. Due to our hypotheses (A3), (A4), we thus lose any control on possible volume preserving turbulences v ∈ V . This, however, has no influence on the system (1.1)-(1.3), which is the subject of our interest here. Inequality (2.24) shows that if U is small in agreement with hypothesis (A1), then also v can be chosen in such a way that hypothesis (A1), interpreted in terms of H 1 , is not violated.

Energy and entropy
In terms of the new variables θ, U, χ , the densities ̺ 0 e, ̺ 0 s of energy and entropy can be written as The energy functional has to be supplemented with the boundary energy term as well as with the gravity potential −̺ 0 gx 3 U . The energy and entropy balance equations now read d dt The entropy balance (3.5) says that the entropy production on the right hand side is nonnegative in agreement with the second principle of thermodynamics. The system is not closed, and the energy supply or the energy loss through the boundary is given by the right hand side of (3.4).
We prescribe the initial conditions for x ∈ Ω , and compute from (3.1)-(3.2) the corresponding initial values e 0 , E 0 Γ , and s 0 for specific energy, boundary energy, and entropy, respectively. Let E 0 and S 0 denote, respectively, E 0 = Ω ̺ 0 e 0 (x) dx, S 0 = Ω ̺ 0 s 0 (x) dx. From the energy end entropy balance equations (3.4), (3.6), we derive the following crucial (formal for the moment) balance equation for the "extended" energy ̺ 0 (e −θ Γ s),θ Γ being a suitable positive constant: We assume that both c(χ) and λ(χ) are bounded from above and from below by positive constants. The growth of s 1 (θ) is dominated by e 1 (θ) as a consequence of the inequality We will use the relation (3.10) to get an upper bound for the solution on the whole time interval (0, ∞). From the identity for all θ, a, b > 0, it follows that we find a constant C > 0 independent of t such that for all t > 0 we have

Main results
We construct the solution of (1.1)-(1.3) by a combined truncation and time discretization scheme. The method of proof is independent of the actual values of the material constants, hence we choose for simplicity We consider the following assumptions on the data.
The liquid phase does not persist for very large temperatures and the behavior of c 1 (θ) as θ → ∞ thus cannot be experimentally verified. We nevertheless believe that the growth condition (ii) in Hypothesis 4.1 is not completely meaningless taking into account the fact that in the interval between 273 and 373 K (0-100 • C), the function c 1 (θ) is convex with a minimum at 35 • C 1 . We introduce the following notation: where (4.5) is to be satisfied for all test functions w ∈ W 1,2 (Ω) and a.e. t > 0, while (4.6)-(4.7) are supposed to hold a.e. in the space-time cylinder that we denote Ω T := Ω × (0, T ) for T > 0, Ω ∞ := Ω × (0, ∞).
In this section we prove the following existence and uniqueness result.

Remark 4.3.
Let us note that we could prove our existence result assuming that κ = κ(θ, χ) = k 1 (θ)k 2 (χ) with the same techniques. Moreover, also uniqueness would hold true in case κ = κ 1 (θ) with an appropriate modification of the boundary condition by means of the standard Kirchhoff transformation technique.

Existence proof
We proceed as follows: We truncate from above the functions depending on θ in (4.5)-(4.7), and discretize the system in time. For the discrete system, we derive upper and lower bounds that enable us to let the time step tend to 0 and prove the existence of a solution to the truncated problem. Finally, we prove a time dependent lower bound and a uniform (in time and w.r.t. the truncation parameters) upper bound on θ , so that the truncation can be removed, and this will conclude the proof of existence of solutions.
It is easy to see then that the latter semi-implicit scheme has a unique solution. Indeed, at each time step, we assume that θ k−1 , U k−1 , χ k−1 are known, and find U k , χ k satisfying (5.13)-(5.14). For τ sufficiently small, (5.13)-(5.14) is an algebraic system for (U k , χ k ) of the form Φ(U k , χ k ) = Y k with Φ : R 2 → R 2 strictly maximal monotone, hence it admits a unique solution. Finally, we insert U k and χ k in (5.12) and solve the resulting coercive elliptic equation, obtaining in that way the desired solution (θ k , U k , χ k ).
Then, we note that the total energy balance still holds true for the discrete system. Indeed, we take (5.12) with w = 1, and denote E k = e R 1 (θ k ) − f 1 (θ c ). We have and using the fact that E k ≥ 0 and that c is convex (cf. Hypo. 4.1 (i)), we get Using now the convexity of λ (cf. Hypo. 4.1 (iii)), we get Hence, from (5.12), using (5.18) and (5.19), we obtain Summing now (5.20) over k = 1, . . . , m , 1 ≤ m ≤ n , we get , with a constant C(T ) independent of τ and R , we check that the left hand side of (5.21) is bounded independently of τ and R . Consequently, all terms in Eq. (5.13) are bounded by a multiple of (1 + B(R)). Similarly, multiplying (5.14) by • χ k and using the fact that • χ k ξ k ≥ 0 for all ξ k ∈ ∂I(χ k ), we obtain the estimates a.e. (5.22)

Lower bound for θ k
Here we derive a lower bound for the approximated absolute temperature θ k . We first rewrite (5.12) for w ∈ W 1,2 (Ω), w ≥ 0 a.e., using (5.13)-(5.14), in the form where we have used again the fact that by definition of the subdifferential. The right hand side of (5.23) is bounded from below by a negative multiple (depending on R ) of θ k−1 θ + k−1 . We can now choose c R in (5.12) sufficiently large in order to get the following inequality for all w ∈ W 1,2 (Ω), w ≥ 0 a.e.: We now compare this inequality with the constant decreasing sequence {v k } defined recurrently as We write (5.25), adding the zero term −div (k(χ k−1 )∇v k ), in the form Assume that θ k−1 ≥ v k−1 (this is true for k = 1). For ε ց 0, (5.28) yields θ k ≥ v k , and by induction we get θ k ≥ v k > v n for all k = 1, . . . , n . By (5.25), we have . This concludes the proof of the lower bound for θ k .

Estimates
Now, we perform the estimates we need in order to pass to the limit as τ ց 0 in (5.12)-(5.14). The right hand side of (5.12) is bounded from above, by virtue of (5.22), by C(T, R)(θ k−1 +1), where C(T, R) is, here and in the sequel, any sufficiently large constant depending only on T and R , and independent of k and τ . Testing (5.12) by w = θ k − θ k−1 , we obtain Using the lower bound for θ k , and choosing χ −1 = χ 0 , we get where C 1 is a positive constant depending on T but not on τ . The elementary inequality θ 2 k − θ 2 k−1 ≤ 1 2τ (θ k − θ k−1 ) 2 + τ 2 (θ k + θ k−1 ) 2 enables us to rewrite the above inequality in the form and C := C(T, R). Inequality (5.29) is equivalent to which yields holding true for τ ≤ 1/(3C). We conclude for all m = 1, . . . , n that Then, we introduce the piecewise constant and piecewise linear interpolants, for t ∈ [(k − 1)τ, kτ ), k = 1, . . . , n , by the formula with a similar notation for U , χ , θ Γ , and P 0 . In particular, we set The previous estimates give immediately that θ (τ ) t bounded in L 2 (0, T ; L 2 (Ω)) , Letting τ tend to 0 and passing to subsequences if necessary, we get the convergenceŝ Now we estimate ∇χ k and ∇U k as follows. From Eq. (5.13) it follows that and, analogously, from (5.14), we obtain Summing up the two previous inequalities, we get We are again in the situation of Eq.
Hence, by (5.31), we obtain Using now the previous estimate on ∇θ τ , we get We already know thatÛ so that the convergences (5.34) take place also for U and χ . We now rewrite (5.12)-(5.14) in terms of the functions θ (τ ) ,θ (τ ) ,ê (τ ) , χ (τ ) ,χ (τ ) ,χ (τ ) ,Ū (τ ) ,Û (τ ) ,θ (τ ) Γ ,P (τ ) 0 . The above estimates allow us to pass to the limit as τ ց 0 and obtain a solution for the following truncated problem The next step consists in proving that θ remains uniformly bounded also from above independently of R , so that the truncation does not become active if R is sufficiently large. The argument is based on the following counterpart of the extended energy balance (3.10), which holds for every solution to (5.35)-(5.37) and every t ∈ (0, T ).

Uniform upper bound for θ
We choose R large enough such that . We may therefore argue as at the end of Section 3 and obtain from (5.39) for all t ∈ (0, T ) that In order to perform the Moser iteration scheme on θ as in [11,Prop. 3.6], we need first to estimate U and U t in terms of θ . Rewriting (5.36) as where, by virtue of (5.40), G(x, t) is bounded above by a positive constant G 0 . Denotinĝ λ(x, t) := t 0 λ(χ(x, s)) ds , we obtain the formula Using Hypo. 4.1 (iii), we get the estimate Now we are ready in order to start the Moser iteration scheme. Choose in (5.35) w(x) = u p , u = ψ R (θ) := (Q R (θ) − R) + , with any p > 1 and with R larger than the constants in Hypothesis 4.1. Then Then, we can rewrite (5.43) as We now prove that the last integral in (5.44) is non-positive if R is sufficiently large. First of all let us note that, if χ t = 0 then it vanishes. Hence, let us consider the case χ t = 0. Then, from (5.37) it follows that χ t (γ(θ)χ t + B R (χ, θ) + C(U, χ)) = 0, hence The last integral in (5.44) is of the form − Ω 1 γ(θ) I 1 × I 2 dx, where We can now estimate from below the last term as follows We have I 2 = 0 if θ ≤ R , while for θ > R we have by Hypothesis 4.1 (i) By virtue of (5.22), we have |C(U, χ)| ≤ C(B 2 (R) + 1). Referring to (5.11), we conclude that there exists R 0 > 1 larger than all constants in Hypothesis 4.1 such that for R ≥ R 0 we have in (5.44) Let us fix now R > R 0 and continue the Moser estimate, rewriting (5.44) as follows We have c(χ)E R P (u) ≥ c * c * p+1 u p+1 , κ(χ) ≥ κ * . Integrating (5.45) in time, we obtain, using Hypo. 4.1 (i),(ii),(iv), as well as the estimates (5.41)-(5.42) and the fact that u(x, 0) ≡ 0, that c * c * p + 1 Ω u p+1 (x, t) dx + 4pκ * (p + 1) 2 The function r(x, t) = CA(U, χ, x, t), where C is a suitable constant, has norm in L ∞ (0, T ; L 2 (Ω)) bounded independently of R by virtue of (5.40). Note that Q R (θ) ≤ u + R . Hence, the function v := u/R satisfies for all p > 1 the inequality The argument of [10,Prop. 4.5] yields v L ∞ (Ω T ) ≤C with a constantC independent of R and T . Consequently, Choosing R sufficiently large such that B(R) > (1 +C)R , we can remove the truncation from (5.35)-(5.37), concluding in this way the proof of existence of a bounded solution to (4.5)-(4.7). If moreover (3.13) holds, then r ∈ L ∞ (0, ∞; L 2 (Ω)), and the upper bound holds globally in Ω ∞ . Indeed, the lower bound for θ in Subsection 5.2 is independent of the time step τ and is preserved when τ ց 0.