Well-posedness for a non-isothermal flow of two viscous incompressible fluids

This work is concerned with a non-isothermal diffuse-interface model which describes the motion of a mixture of two viscous incompressible fluids. The model consists of modified Navier-Stokes equations coupled with a phase-field equation given by a convective Allen-Cahn equation, and energy transport equation for the temperature. This model admits a dissipative energy inequality. It is investigated the well-posedness of the problem in the two and three dimensional cases without any restriction on the size of the initial data. Moreover, regular and singular potentials for the phase-field equation are considered.


1.
Introduction. The interfacial dynamics of a mixture of different fluids plays an important role in the hydrodynamic theory due to the increasing science and engineering applications. In the classical approach, the interface is assumed to have zero thickness and this leads to a free boundary problem. Even though classical models have been applied successfully to many situations, they break down, for instance, when parts of the interface develop singularities. Moreover, in many real processes the interfaces are actually not sharp. Diffuse-interface models emerge as an alternative description of the interface to handle these plights. The main idea behind this approach is to take into consideration from the beginning that physical interfaces have some thickness, maybe small, and also a structure. See [2] for a review on diffuse-interface models.
This work is concerned with a non-isothermal diffuse-interface model which describes the motion of a mixture of two viscous incompressible fluids. The fluids are assumed to have matched densities and the same viscosity and thermal conductivity. The model consists of modified Navier-Stokes equations coupled with a phase-field equation given by a convective Allen-Cahn equation, and a equation for the temperature. More precisely, we consider the following system u t + u · ∇u − ∇ · (ν(θ)Du) + ∇p = (− ∆φ + F (φ))∇φ − α∆θ∇θ (1) θ t + u · ∇θ = k∆θ (4) in Ω × (0, ∞), where Ω is a bounded domain in R n , n = 2, 3, with smooth boundary ∂Ω. Here, u denotes the mean velocity of the fluid mixture, p is the pressure, the phase-field variable φ represents the volume fraction of the two components, θ is the temperature, ν(·) ≥ ν 0 > 0 is the viscosity of the mixture, > 0 is a small parameter related to the interfacial thickness, F (φ) is the potential energy density, α > 0 is associated to the thermal expansion coefficient, γ > 0 is the relaxation time of the interface, and k > 0 the thermal conductivity. Du = 1 2 (∇u + (∇u) T ) corresponds to the symmetric part of the velocity gradient.
For the potential energy density, we treat both cases regular and singular. It includes the thermodynamically relevant potential of logarithmic type where κ 1 , κ 2 > 0, and the usual double-well potential which is often employed as a polynomial approximation for the logarithmic potential, [5,3].
We supplement system (1)-(4) with initial and boundary conditions u(0) = u 0 , φ(0) = φ 0 , θ(0) = θ 0 in Ω, where η is the outward normal on ∂Ω. The aim of this paper is to prove the well-posedness to problem (1)- (6) in both cases, when F is a regular and a singular potential. The precise definition of the potential will be introduced in Section 2.
In the isothermal case, the model (1)-(4) reduces to the Navier-Stokes-Allen-Cahn system investigated in [19,20,39]. A simplified model, by replacing the Navier-Stokes equations by the Stokes equations, was analysed in [14], where numerical simulations show that the system captures the basic features of a binary fluid behaviour. Several models that include an Allen-Cahn type equation for the phase-field have been considered in the literature. We can mention: [27] for the two fluids problem with surface tension, [9] for vesicle dynamics, [40] for drop formation processes, [4,29] for binary alloys, [25,7] for nematic liquid crystal problems, [12] for compressible two viscous fluids, [24,23] for two fluids with non-matched densities, among others.
A further diffuse-interface model for binary fluids, which consists in coupling the Navier-Stokes equations with the Cahn-Hilliard equation, was investigated by several authors, see, for instance, [1,18]. As pointed out in [26], from the energetic point of view, the choice between the Cahn-Hilliard equation and the Allen-Cahn equation determines the special dissipative mechanism of the system. From the mathematical point of view, the analysis of the system based on the Allen-Cahn equation is more delicate due to the lower regularity of the phase-field; however, in this case, a maximum principle is available which does not hold for the Cahn-Hilliard equation.
On the other hand, it is not trivial question to include temperature dependence in a way such that the obtained models are at the same time thermodynamically consistent and mathematically tractable. In the case of a double-well potential, [33] introduced a non-isothermal model for a mixture of two fluids with thermoinduced Marangoni effects, for which the mathematical analysis was performed in [37,38]. In such models, the dissipative energy inequality for the system is valid under a smallness condition on the initial temperature. Lets us mention very recent works [10,11] where the authors proposed and investigated a non-isothermal Navier-Stokes-Cahn-Hilliard system. That model contains a highly non-linear temperature equation for which they are able to obtain only weak/entropy solutions.
One of the main advantages of the model under consideration is that, for any non-negative potential, it satisfies a dissipative energy inequality for any initial data, which it is expected from the physical point of view. Moreover, in the general case, the energy estimate allows to obtain the existence of global solutions. Furthermore, comparing with previous non-isothermal models, our model has higher order nonlinearities in the fluid equations which bring several new difficulties to an already hard problem.
It is important to mention that the previous papers on non-isothermal two-fluids problem [37,38,10,11] considered only the double-well potential. The analysis for singular potential requires the phase-field to belong to the physically relevant interval (−1, 1). The main difficulty, when considering such singular potentials, is to ensure that the phase-field remains strictly separated for all time from the singular values of the potential. For related results concerning singular potentials, we can cite [6] for Caginalp phase-field system, [19] for a two-dimensional model for an isothermal incompressible two-phase fluids, and [12] for a phase-field model for two-phase compressible fluids.
The structure of the paper is as follows. In Section 2, we introduce the notation and state the main results of this paper. In Section 3, we show the existence of weak solutions for our system with F being a non singular potential for any initial data, by using the semi-Galerkin method. The maximum principle for the phase-field and temperature equations are important to overcome the difficulties that come from the higher order non-linearities. In Section 4, we prove the existence and uniqueness of the strong solution with F a non singular potential. We show the existence of global strong solution for n = 2 and local for n = 3. In Section 5, we treat the singular potential case. Finally, in Section 6 we briefly comment on some modeling aspects of the problem.
2. Main results. Let us firstly introduce some notation which will be used all through the paper. Let Ω be a bounded domain in R n , n = 2, 3, with smooth boundary ∂Ω. We denote by |Ω| the measure of this set. For 1 ≤ p ≤ ∞, L p = L p (Ω) stands for the Lebesgue space and H m = W m,2 (Ω), 0 ≤ m < +∞, for the Sobolev space.
For the sake of simplicity, · will denote the L 2 or (L 2 ) n norms and (·, ·) the inner product in R n and the tensor product in R n 2 ; the norms on H m and (H m ) n are indicated by · H m ; we write u ∈ L 2 , even when u is a vector field, meaning that all of its components are in L 2 , and so on. The difference will be clear from the context.
We introduce the standard functional spaces of divergence-free vector fields Observe that, by Poincaré and Korn inequalities, ∇v and Dv are equivalent norms in V . The duality product between V and V will be denoted by ·, · . Let P be the orthogonal projection from L 2 onto H. We shall denote by {w i } and {λ i }, respectively, the eigenfunctions and eigenvalues of the Stokes operator A = −P ∆ defined in V ∩ H 2 . It is well known that {w i } are orthogonal and complete in the spaces H, V and V ∩ H 2 (see [36]). Notice that, from the Helmholtz decomposition of −∆v, we can write −∆v = Av + ∇q for v ∈ H 2 ∩ V and some q ∈ H 1 with Ω q = 0. Moreover, from Lemma 3.4 in [28], there exist C > 0 and, for any δ > 0, C δ > 0 such that We recall here Ladyzhenskaya and Gagliardo-Nirenberg interpolation inequalities that will be frequently used in the text (see [36,17,31]). For any v ∈ H 1 that vanishes on the boundary of Ω it holds and for u ∈ L q ∩ W m,r , 1 ≤ q, r ≤ ∞, for some positive constant C 1 with the only exception: if 1 < r < ∞ and m − j − n r is a non-negative integer, then the above inequality holds only for j m ≤ α < 1. Throughout this paper, C will denote a positive constant which may vary from line to line.
Next, we state our assumptions on the potential energy density. Motivated by the double-well potential, F will be called non singular when F ∈ C 2 (R) and satisfies the following conditions F (1) ≥ 0 and F (−1) ≤ 0.
Note that the classical double-well potential F (s) = 1 4 (s 2 − 1) 2 satisfies (11)- (12). Let us also observe that assumption (12) guarantees the maximum principle for the phase-field equation. Moreover, it is not difficult to see that condition (11) implies that there exist some positive constants c i , i = 1, 2, 3 such that where F (s) = s 0 F (r)dr.
Observe that with above conditions, property (13) is still true for all s ∈ (−1, 1). We introduce the concept of weak solution to problem (1)-(6).
Additionally, we say that (u, φ, θ) is a global weak solution if its restriction to (0, T ) is a weak solution for any T > 0.
We observe that, in the weak formulation, since the vector field v is divergence free and v = 0 on ∂Ω, the terms involving ∇p and F (φ)∇φ = ∇F (φ) vanish.
Next, we introduce the notion of strong solution to problem (1)-(6).
We now state our main result about existence of global weak solutions for regular potentials.
Then problem (1)-(6) has at least one global weak solution that satisfies . Ω × (0, ∞), and the following energy inequality, a.e. t > 0, With stronger assumptions on the data, we can prove the existence and uniqueness of strong solutions for regular potentials.
Then problem (1)-(6) has a unique strong solution which is global for n = 2 and local for n = 3. Moreover, it depends continuously on the initial data in H × H 1 × H 1 .
Let us notice that the uniqueness of the strong solution will be a consequence of the continuous dependence of the solution with respect to the initial data.
Finally we can extend the previous theorems to singular potentials.
Theorem 2.5. Suppose that the singular potential F satisfies (14). Theorems 2.3 and 2.4 remain true if we assume that initial condition for the phase-field satisfies φ 0 L ∞ < 1.
We observe that the size of the constants , α, γ and k do not play any role in the proof of the main results. Therefore, we will assume, without any loss of generality, that = α = γ = k = 1.
3. Weak solutions. In this section we prove Theorem 2.3 on the existence of weak solutions to problem (1)-(6) for regular potentials. To this end, we apply the semi-Galerkin method. Let w i and λ i , i ∈ N, be the eigenfunctions and eigenvalues of the Stokes operator A. For all m ∈ N, we denote by P m the orthogonal projection from H onto H m = span{w 1 , . . . , w m }.
Let us fix T > 0. For every m ∈ N we consider the following approximate problem of finding The existence and uniqueness of the local solution defined on (0, T m ) for some T m > 0 can be obtained by a fixed point argument. The proof is done in Proposition A.3 of the Appendix. We note that φ m and θ m satisfy the maximum principle, more precisely, |φ m | ≤ 1 and |θ m | ≤ θ 0 L ∞ a.e. in Ω × (0, T m ). Moreover, by taking u m as test function in (16), multiplying equation (17) by −∆φ m + F (φ m ), equation (18) by −∆θ m , integrating in Ω and adding the resulting identities, we find that where we have used that (u m · ∇u m , u m ) = 0 and (u m · ∇φ m , F (φ m )) = 0. Hence, integrating in time and using that Thus, we can take T m = T . Moreover, we infer from the above estimates, after integration in time and using standard elliptic estimates, (see [21, independently of m. From the uniform estimates we deduce that there exist functions u ∈ L ∞ (0, T ; H) In order to pass to the limit on the non-linear terms we need some strong convergences. To this end, we will estimate the derivative with respect to t of the approximate solution.
Notice that by the identity ∇ · (∇ω ⊗ ∇ω) Therefore, from (16) and using the maximum principle for θ m , we have that To estimate the right hand side we use the Gagliardo-Nirenberg interpolation inequality (10) and the fact that {φ m } m and {θ m } m are bounded in L ∞ (0, T ; L ∞ ). In particular, we can bound and, similarly, By using the Ladyzhenskaya inequality (9), when n = 3 we have that From equation (17) there follows Hence we can pass to limit in (16) , and find out that there exists at least one weak solution for (1)-(6) in (0, T ).
Lets show the energy inequality (15) which allows to extend the solution to (0, ∞). By integrating (21) in time, it follows that Since the right hand side does not depend on t, we have that the approximate solution is defined in (0, ∞) and taking the limit as m goes to ∞, we see that (15) holds.

Remark 1.
It is worth to mention that, due to elliptic estimates for the Neumann problem (see, for instance, [21]) and the estimates of φ and θ in L 2 loc (0, ∞; L 2 ) there follows that φ and θ belong to L 2 loc (0, ∞; H 2 ).
4. Strong solutions. In this section we prove Theorem 2.4 about existence and uniqueness of strong solutions to problem (1)- (6). In order to do this, we will obtain a differential inequality for higher order norms for the approximate solution.
To simplify the notation, we will omit the superscript m in the following. First, we take Au ∈ H m , where A is the Stokes operator, as test function in (16) to find 1 2 Recalling that −∆u = Au + ∇q we can write and using the fact that Au = 0 on ∂Ω and ∇ · Au = 0, there follows Therefore, Next, we apply the gradient to equation (17) for φ and to equation (18) for θ, then multiply by −∇∆φ and −∇∆θ, respectively, and integrate over Ω. Since ∂φt ∂η = ∂θt ∂η = 0 on ∂Ω, we have that −(∇φ t , ∇∆φ) = 1 2 d dt ∆φ 2 and analogously for θ. Finally, we add the resulting identities to the previous one to obtain We will estimate the right hand side of (23) term by term. These bounds will depend on the dimension and yield a different differential inequality for n = 2 and n = 3. Thus, in the next Subsection, we treat the two-dimensional case and we are led to the existence of global solutions. In Subsection 4.2 we deal with the three dimensional case and show the existence of local solutions. In the last Subsection we prove the continuous dependence of the strong solution with respect to the initial data.
4.1. Global solutions for n = 2. The proof is divided in two steps. In the first step we derive a nonlinear differential inequality which allows to show the existence of local strong solutions. The main difficulty arises from the term J 5 involving the function q. Let us notice that if ν is constant, this term will no appear and we could obtain directly the existence of a global strong solution. In the second step, by using the regularity of the Stokes system with variable coefficients, we obtain another differential inequality which allows to extend the local solution.
Then, we have, for all t ≤ T m , Integrating to 0 to t, we obtain for all t ∈ [0, T m ], which implies that G(T m ) + 1 < 2µ. Thus, we arrive at a contradiction. Thenceforth, T m = T * and, by elliptic estimates, we conclude It is not difficult to see that This shows the existence of local strong solutions to (1)-(6).
As in (28), we estimate q by using Gagliardo-Nirenberg inequality (10) and (7); thus, for any δ > 0 we have that where C is independent of ε and δ.
Finally, the bound (29) for J 8 does not depend on the dimension.
The existence of a local strong solution follows as in the the two-dimensional case.

Continuous dependence and uniqueness.
In this subsection we show the continuous dependence of the strong solution with respect to the initial data which yields the uniqueness of the solution.
We will derive a differential inequality for some norms of u, φ and θ which allows to apply the Gronwall Lemma in order to conclude the result.
Observe that u, φ and θ satisfy θ t + u · ∇θ 1 + u 2 · ∇θ = ∆θ a.e. in Ω × (0, T ), and in Ω. By taking u as test function in (38), the scalar product in L 2 of equation (39) with φ and with −∆φ, and of equation (40) with θ and with −∆θ, respectively, adding up the resulting identities and arranging terms, we arrive at 1 2 For the sake of simplicity, we proceed to estimate the right hand side of (41) term by term in the three-dimensional case. Similar result can be obtained in the two dimensional case with minor modifications.
By using Hölder, Ladyzhenskaya (9) and Young inequalities, we have where ε > 0 is a constant that will be determined later. By using Gagliardo-Nirenberg inequality (10) and elliptic estimates we can bound Thus, as before, by using Hölder and Young inequalities,

Analogously, there follows
The next two terms can be estimated by using the Sobolev embedding H 1 ⊂ L 4 and the Poincaré inequality, For I 6 and I 7 , we estimate ∇φ L 4 using (42). Hence, For the last three terms, we use the Mean Value Theorem together with the fact that φ and θ belong to L ∞ (0, T ; L ∞ ) to estimate Thus, by plugging the above estimates in (41) and taking ε > 0 small enough, we have that The regularity of the strong solution implies that I(t) is integrable in [0, T ], thus Gronwall Lemma gives the continuity dependence of the solution with respect to the initial data in H × H 1 × H 1 . The uniqueness follows in a straightforward manner.
This completes the proof of Theorem 2.4.

Remark 2.
We observe that in the two dimensional case it is possible to show the uniqueness of weak solutions when the viscosity coefficient is constant. To this end, it only necessary to re-estimate I 1 , I 2 , I 3 , I 6 and I 7 in (41). This can be done in an analogous way but now using Ladyzhenskaya inequality (8) and Gagliardo-Nirenberg inequality for n = 2.
5. Singular potential case. All the results on existence of weak and strong solutions to problem (1)-(6) can be extended to singular potentials. Indeed, it is enough to show that the phase-field does not coincide with −1 and 1, the so-called pure phases, thus F will remain strictly separated from its singularities. Here, we assume that φ 0 L ∞ < 1. We recall that for singular potentials, F ∈ C 2 (−1, 1) satisfies lim s→±1 F (s) = ±∞ and lim s→±1 F (s) = +∞.
We first state the maximum principle for the phase-field φ when the potential is singular that can be proved as Theorem 6.1 in [19]. In fact, the proof uses the fact that if φ 0 L ∞ < 1 and F satisfies (43), then there exists δ ∈ (0, 1) such that Whence it is deduced that |φ(x, t)| ≤ δ. So, we have then φ is strictly separated from the values −1 and 1, that is, there exists δ ∈ (0, 1), depending only on φ 0 L ∞ , such that Next, we introduce a new potential that is an approximation of F . Motivated by [6,19], we define where δ is the constant given in (44).
Let F δ (s) = s 0 F δ (r) dr and observe that properties (13) are true for all s ∈ R: where c i > 0, i = 6, 7, 8, are independent of δ. Thus, F δ satisfies the properties of a regular potential. Now, we consider problem (1)-(6) with F δ instead of F and construct an approximate problem for this new system by using the semi-Galerkin method. The proof of existence and uniqueness of the solution (u m δ , φ m δ , θ m δ ) to the approximate problem can be done by a fixed point argument as in Proposition A.3 of Appendix.
Let us notice that the approximate equation of φ m δ has a unique strong solution and satisfies the maximum principle (Proposition 5.1). Hence, φ m δ (t) L ∞ ≤ δ < 1, for all m. Thus, we have F δ (φ m δ ) = F (φ m δ ), so the approximate problem with F δ coincides with the approximate problem with F , the singular potential. Therefore, we conclude that (u m δ , φ m δ , θ m δ ) = (u m , φ m , θ m ) is the solution for the approximate problem of the original system (1)- (6).
Since φ m (t) L ∞ ≤ δ < 1 for all t ≥ 0 and for all m, the proof of existence of weak and strong solutions follows as in the regular potential case. 6. Interpretation of the model. In this section, we comment some modeling aspects of the problem which show that the model admits a dissipative energy law (similar arguments can be found, for instance, in [26,40]).
Consider a physical domain Ω filled with two incompressible fluids taking into account small fluctuations of the density due to thermal expansion (Boussinesq type approximation) for which the Lyapunov functional or "free energy" can be expressed in the form where φ is the phase-field, θ the temperature, 2 |∇φ| 2 corresponds to interfacial accumulated energy (non-local interactions between the components), > 0 is a small parameter related to the interfacial thickness, F (φ) is the non singular or singular potential (the bulk part of the mixing energy), α 2 |∇θ| 2 is the energy due to the variation of the specific volume (inverse of the density) coming from thermal expansion, and α > 0 is a parameter associated to the thermal expansion coefficient.
We notice that gradient of temperature effects in the free energy have been considered in a different context, see, for instance, [15,22].
According to the energetic variational approach, as in [26,13,40], the evolution equation for the phase-field φ is then given by where D t = ∂ t + u · ∇ is the convective derivative being u is the mean velocity of the fluid mixture, γ > 0 the elastic relaxation time and δE δφ represents the variation of the energy with respect to φ. So, we arrive at a convective Allen-Cahn type equation (45) The balance of the linear momentum for homogeneous incompressible Newtonian viscous fluids, assuming no external forces, has the following differential form where ρ is the density, σ is the stress tensor, and f micro is an internal force due to interfacial interaction and thermal fluctuations. The viscous stress tensor is given by where p is the hydrostatic pressure and ν(θ) ≥ ν 0 > 0 is the viscosity coefficient which can depend on the temperature. Hence, assuming that the density ρ = 1 for simplicity and using the fact that ∇ · u = 0, we find The Least Action Principle (Principle of Virtual Work [8,16]) can give an expression for f micro for some function q. By absorbing the pure gradient term into the pressure we arrive at the following modified Navier-Stokes equation for the velocity field ∇ · u = 0.
Finally, the balance of energy can be written as where κ > 0 is the thermal conductivity. Hence, we find the internal energy equation describing the evolution of temperature where k = κα > 0. We close the system (45)-(47) by adding the initial conditions, no-slip boundary condition for the velocity u and no-flux boundary conditions for the phase-field and temperature, that is, where η is the outward normal on ∂Ω.
Considering a Lyapunov functional (the total energy of the system) E tot as the sum of the kinetic energy of the fluid and the interfacial energy E(φ, ∇φ, ∇θ), that is, it is not difficult to see that, for smooth enough solutions, it holds Thus, the system admits a dissipative energy law for any non-negative potential F (φ). For general potentials, this energy allows to show the existence of global solutions.
Appendix. In this section we prove the existence and uniqueness of the local solution to the approximated problem (16)- (20). First we recall some facts concerning the phase-field and temperature equations with smooth coefficients.
Given v m ∈ C([0, T ]; H m ), consider following problem The existence and uniqueness of solution to this problem is obtained by a standard Galerkin method. Moreover, the maximum principle is true (see, for instance, [25,38,7]). We resume these results in the following Lemma. and φ m t ∈ L ∞ (0, T ; L 2 ) ∩ L 2 (0, T ; H 1 ). Similarly, for the temperature equation we have the following result where M is a large enough constant that will be chosen latter. Now, given v m ∈ C([0, T ]; H m ) as above, by Lemma A.1, there exists a unique solution φ m to problem (48) which satisfies |φ m (x, t)| ≤ 1 a.e. in Ω × (0, T ). Moreover, by multiplying equation (48) Similarly, by Lemma A.2, there exists a unique solution θ m to problem (49) which satisfies |θ m (x, t)| ≤ θ 0 L ∞ a.e. in Ω × (0, T ) and the following estimate ∇θ m 2 ≤ e CmM T ∇θ 0 2 .
From the local existence theorem of the classical ODE theory, there exists a unique solution ( g m k (t)) 1≤k≤m defined on [0, T m ), for some T m > 0 depending on m, and g m k ∈ H 1 (0, T m ). We will prove that T m = T . To this end, by taking u m as test function in (16) To estimate the right hand side, we use the identity ∇ · (∇ω ⊗ ∇ω) = ∇ |∇ω| 2 2 + ∆ω∇ω, where ⊗ is the Kronecker product, which is the matrix with entries (i, j) are given by (a ⊗ b) ij = a i b j , to write −(∆φ m ∇φ m , u m ) = (∇φ m ⊗ ∇φ m , ∇u m ), −(∆θ m ∇θ m , u m ) = (∇θ m ⊗ ∇θ m , ∇u m ).
As H m is a finite dimensional space, it holds ∇u m L ∞ ≤ C m ∇u m for a constant C m depending on m. Thus, we obtain the following estimates Plugging the previous estimates in (53) and using estimates (51)-(52) of φ m and θ m , we arrive at d dt u m 2 + Du m 2 ≤ C m .
Integrating from 0 to t ≤ T m , we find Therefore, we conclude that T m = T and that u m is bounded in L ∞ (0, T ; H) ∩ L 2 (0, T ; V ) by a constant that depends on m and M . From the EDO system, we see that d dt g m i is bounded in L 2 (0, T ) and thus u m is bounded in H 1 (0, T ; H m ). To show that v m = u m we apply the Schauder Fixed Point Theorem. Lets define the following operator The proof of the uniqueness of the solution is standard, so we omit the details. The further regularity for φ m and θ m follows from Lemmas A.1 and A.2. This completes the proof.