TWO-PHASE INCOMPRESSIBLE FLOWS WITH VARIABLE DENSITY: AN ENERGETIC VARIATIONAL APPROACH

. In this paper, we study a diﬀuse-interface model for two-phase incompressible ﬂows with diﬀerent densities. First, we present a derivation of the model using an energetic variational approach. Our model allows large density ratio between the two phases and moreover, it is thermodynamically consistent and admits a dissipative energy law. Under suitable assumptions on the average density function, we establish the global existence of a weak solution in the 3D case as well as the global well-posedness of strong solutions in the 2D case to an initial-boundary problem for the resulting Allen-Cahn-Navier-Stokes system. Furthermore, we investigate the longtime behavior of the 2D strong solutions. In particular, we obtain existence of a maximal compact attractor and prove that the solution will converge to an equilibrium as time goes to inﬁnity.


1.
Introduction. This paper is devoted to the study of the following Allen-Cahn-Navier-Stokes (ACNS) system which models hydrodynamics of the two-phase incompressible flows with different densities: Here, φ is the phase field function that labels different species, u is the velocity and p is the pressure. µ, λ and γ are positive physical parameters. ρ(·) is the average density which is a given function of φ. The physical relevant energy density function f that represents the two phases of the mixture usually has a double-well structure. Without loss of generality, we assume that with some small parameter ε > 0. Let Ω ⊂ R n (n = 2, 3) be a bounded domain with smooth boundary. We complement the above system with the following boundary conditions: u ∂Ω = 0, ∂ ν φ ∂Ω = 0 (5) and initial data: The study of the interfacial dynamics of immiscible and incompressible twophase flows is of great importance in the hydrodynamic theory of complex fluids [13,16]. Recently, the diffuse-interface approach, or sometimes called the phase field method, becomes a very successful and major tool to study a variety of interfacial phenomena [2,3,4,12,15,25,31,26,38,32]. In this framework, a labeling function is introduced to identify different species while the sharp interfaces are thus replaced by narrow transition layers, whose thickness is approximately characterized by the small parameter ε > 0. Within this "thin" transition region, the fluid is mixed and stores certain amount of "mixing energy" (see [8]). Instead of tracking the interfaces explicitly, the dynamics of interfaces (now recognized as zero level sets of the order parameter) can be simulated on a fixed grid and the topological changes of the interfaces can be described in a natural way [27].
In recent years, many researchers have employed the phase field approach in various fluid environments [2,3,4,15,25,31,26,38] and have carried out extensive analysis and numerical studies on the two-phase flows in the literature. For phase field model of incompressible binary fluids with identical densities, or with small density ratios where Boussinesq approximation can be applied in practice, we refer the readers to [1,18,19,20,11,39,38] and references therein for detailed derivations and mathematical analysis. However, in most cases, the density differences of the components are not negligible, whence studies on the incompressible two-phase fluids with non-matched densities become even more interesting and challenging. In [27], the authors derived a quasi-incompressible phase model with a mass averaged velocity that resulted in a non-divergence free velocity field, see also [2,32]. On the other hand, alternative phase-field models were derived in [12] and [14] using volume averaged velocity that leads to a solenoidal velocity field. However, it is not clear to us whether these models satisfy any dissipative energy law. It is worth noting that in [3,4], the authors derived two different diffuse interface models for incompressible two-phase flows with different densities, both of which were thermodynamically consistent. Another thermodynamically consistent phase field model of Cahn-Hilliard type for incompressible two-phase flows was also proposed in [26] via energetic variational approach and was studied from a numerical point of view. One objective of the present paper is to derive an Allen-Cahn type model in the same spirit of [26]. To our knowledge, theoretical results on models for incompressible binary flows with different densities are quite limited. In [12], Boyer considered a quasi-incompressible diffuse interface model for fluids with non-matched densities and existence of local strong solutions as well as existence of global weak solution provided that the density difference was sufficiently small. In [5], Abels et al. proved the existence of weak solutions for the model derived in [4] with a singular potential and later in [6], they extended the existence results to the same model with further assuming that mobility was degenerate. Now, let us briefly report the main theoretical results of this paper, where the notation of the spaces is given in Section 3. The first one is about the global solvability of our problem. Theorem 1.1. Let assumptions (23)- (25) hold and T > 0 be a given constant.
in Ω, problem In addition, the strong solution (u(t), φ(t)) will converge to a single equilibrium (0, ψ) in V × H 2 (Ω) with ψ being a solution to the corresponding stationary Allen-Cahn equation (78). In particular, there exists C > 0 and θ ∈ (0, 1 2 ) depending on ψ such that for all t ≥ 0, Remark 1. Our convergence results generalize those in [19] for incompressible binary flow with matching density and we note that the corresponding stationary equations are the same since u decays to zero and the quadratic nonlinear coupling term in (3) vanishes.
Before concluding the introduction part, we would like to stress some new features of the present paper.
First, we propose a novel diffuse-interface model of Allen-Cahn type for the twophase incompressible flow with non-matched densities and a solenoidal velocity field. The way we employed in the present paper is the so-called energetic variational approach (cf., [4,12,14]). In such a framework, the overall dynamics of an isothermal hydrodynamical system are governed by the First Law of Thermodynamics and the Second Law of Thermodynamics that lead to where E total is the sum of kinetic energy and the free energy, is the entropy production (energy dissipation in this case). Therefore, the expressions of the total energy and dissipation functionals and the kinematics relation of the variables should contain all the physical ingredients, assumptions and limitations of the system. Observing the multi-scale nature of our problem, we apply the Least Action Principle (LAP) and the Maximum Dissipation Principle (MDP) to derive the expressions of conservative forces and dissipative forces of the system in both macroscopic scales and microscopic scales. Then the resulting system is followed merely by adding all terms together due to the macroscopic force balance as well as microscopic force balance.
The other difficulty in modeling is dealing with non-matching densities. As the density is a macroscopic quantity, it may be different from the direct average of the microscopic descriptions, such as of the phase field parameters. For example, since we have to take into account the interaction between two different particles, the mixture of two incompressible fluids may not be incompressible any more. Traditionally, there are two ways of modeling for this problem, incompressible and quasi-incompressible. In the former approach, the volume averaged velocity is assumed to be incompressible everywhere [12,14,3,4] while in the latter one, the mass averaged velocity is assumed to satisfy the mass conservation instead of incompressibility [27]. In contrast, one fundamental assumption of our framework is that the incompressible binary fluid is moving at a unified effective velocity both macroscopically and microscopically and the velocity field satisfies incompressibility everywhere [26,25,31]. Then by the energetic variation approach, we derive a highly nonlinear PDE system consisting of an incompressible Navier-Stokes equation with varying density function and a convective Allen-Cahn type equation with presence of a quadratic term of the velocity, i.e., ρ (φ) 2 |u| 2 in (3). This additional nonlinear coupling term represents certain density force that well captures the macroscopic fluid effect on the microscopic descriptions due to density differences. With matched density, this term vanishes and our model recovers the usual phase-field model for binary incompressible fluids.
Besides the modeling, we also present several analytical results for the resulting ACNS system on existence of solutions and their longtime behavior. As we mentioned above, analytical and numerical results for two-phase incompressible flow models that are valid for large density ratio between different species seem quite limited before the present work (cf. [6,26,31]). Most studies of the phase-field models for binary fluids have been restricted to the cases with the same density or with small density differences. In the latter case, Boussinesq approximation can be used, where the variable density is replaced by a background constant density value and an external gravitational force is added to model the effect of density force (see, [25,38]). However, in our model, the density ratio between two phases can be quite large and hence the Boussinesq approximation is no more physically valid. Under the circumstances, a nonlinear term ρ (φ) 2 |u| 2 presents in the Allen-Cahn equation and brings new difficulties in the analysis which also appears in a Cahn-Hilliard type model derived in [3]. For example in 3D, the regularity of φ and φ t is much weaker than in the usual matched density case. Moreover, the density function no longer satisfies the conservation of mass since it is a given function of φ now. As a result, some useful relations between the dependent variables break down in our case and hence the existence theory of density-dependent incompressible Navier-Stokes equations cannot be applied directly (see e.g., [22,35,9]). We need to carry out more delicate estimates and careful compactness argument to prove the existence of 3D weak solutions. Furthermore, instead of employing a singular potential as discussed in [19,5,6], we propose an exterior convexity assumption on the density function outside region [−1, 1] (see (25)), which restricts the value of the phase field parameter in the physically reasonable region [−1, 1] by maximum principle (see Lemma 3.2) and this property plays a crucial role in the whole proof. Last, to our knowledge, longtime behavior of solutions to incompressible two-phase flows with large density ratio is studied for the first time. Based on a dissipative energy law of our system, we are able to prove existence of a maximal attractor and the convergence of global solution towards an equilibrium as times goes to infinity via the so-called Lojasiewicz-Simon approach. Our results recover those in [19] if there is no density difference between the species.
The rest part of this paper is organized as follows. In Section 2, the model is derived via the energetic variational approach. First, we introduce the relevant energies. Then applying the Least Action Principle and the Maximum Dissipation Principle, we obtain a coupled ACNS system for the two-phase incompressible flow. In Section 3, we study the existence of a weak solution to the initial-boundary value problem (1)-(6) based on the Galerkin method and in Section 4 we further prove the well-posedness of strong solutions of our problem. In the last section, we analyze the longtime behavior of the strong solutions as time goes to infinity. More precisely, existence of maximal compact attractors as well as convergence of solutions towards steady states were proved. In the Appendix, we give a proof of an auxiliary lemma for a convective Allen-Cahn equation.

2.
The derivation of the model. In the same spirit of [26] , we derive a diffuseinterface model of Allen-Cahn type for the incompressible binary flow based on the so-called energetic variational approach. Let us consider a mixture of two immiscible, incompressible fluids in a bounded domain Ω ⊂ R 3 with densities ρ 1 and ρ 2 . A "phase" function φ(x, t) is introduced to identify the two fluids such that φ(x, t) = 1 fluid 1, −1 fluid 2, with a thin, smooth transition region of width O( ). Now, let us introduce the basic variable in the context of hydrodynamics, the flow map (particle trajectory) x(X, t), where X is the original labeling (Lagrangian coordinate) of the particle and x is the current (Eulerian) coordinate. For a given velocity field u(x, t), the flow map is defined by The deformation tensor F (X, t) is defined as Then the incompressibility is represented as Notice that when expressed in Lagrangian coordinate φ is independent of time, φ(x(X, t), t) = φ 0 (X), soφ = 0 or, equivalently, in an Eulerian frame. As we mentioned in the previous section, since we consider the hydrodynamics of the mixture as an isothermal system, it follows from the First Law of Thermodynamics and the Second Law of Thermodynamics that the overall dynamics is governed by relation (7). In order to characterize the macroscopic flow properties which self-consistently take into consideration the microscopic properties of the mixture, we define the total energy of our system as follows: Here the first part is the macroscopic kinetic energy and the second part is the Helmholtz free energy consisting of a hydrophilic term describing the tendency of mixing of the fluids and a hydrophobic term accounting for the tendency of separation of the fluids. The constant λ represents the competition between the kinetic energy and the free energy. The dissipation of the energy on the right hand side of (7) is given as where the first part accounts for the macroscopic dissipation due to viscosity and the second part comes from microscopic internal damping during the mixing. Summing up, the energy law of our system now reads: In order to derive the force balance equations from the above energy law, we employ the Least Action Principle for the Hamilton part of the system and the Maximum Dissipation Principle for the dissipative part.
As we mentioned before, our problem is in essence a multi-scale system. Therefore, we need to establish the macroscopic force balance as well as microscopical (potential) balance. First, we derive the macroscopic force balance equation. To this end, we introduce the action function, after Legendre transformation [10], in view of the kinematic transport relation (9) and volume preserving constraint (8) as following: The LAP states that the conservative force is given by the taking variation (derivative) of the action function with respect to displacement. Since our flow map is volume preserving, we use volume preserving diffeomorphisms to perform the variation. Suppose there is a one-parameter family of such flow maps x ε (X, t) satisfying Then we deduce that whereỹ(x(X, t), t) := y(X, t). Now, we can calculate the variation of the action function as follows.
Thus the revertible part of macroscopic force takes the following form: where p 1 is the Lagrangian multiplier because of the volume-preserving constraint and ∇φ ⊗ ∇φ is the stress tensor accounting for the microscopic capillary force due to surface tension at the mixing surface.
On the other hand, MDP, by Onsager and Rayleigh [28,29,30], yield the dissipative force is given by Here the factor 1 2 is consistent with the choice of quadratic form of the "rates" which in turn reflects the linear response theory for long-time near equilibrium dynamics [23,30]. Thus, by taking variation of the dissipation term, we obtain the macroscopic dissipative part of the system: where p 2 is the Lagrangian multiplier for the incompressibility constraint ∇ · u = 0. As a result of the macroscopic force balance, we obtain the momentum equation for the system: where p = p 1 − p 2 . The dynamics of φ is described by the microscopic force balance which can be viewed as a relaxation of the pure transport equation (9). We may apply LAP and MDP with respect to φ andφ, respectively to derive the following convective Allen-Cahn equations (cf. [26]): The right-hand side of the above equation does not satisfy the conventional macroscopic frame-invariant principle (see [7]). In fact, the Allen-Cahn type of dynamics shall not be viewed strictly as a physical law. Rather it is just a relaxation of the pure transport physics. The right hand side is not necessarily a frame-indifference quantity, unlike that in the model proposed by Abels et al. [4] which is the chemical potential. Moreover, it is noticed that our model is in essence a system of macromicro model for which there is no point just to keep the frame-indifference of the macroscopic motions. Although in a macroscopic continuum description, it tries to incorporate the effects from the microscopic structures of the mixture. It is clear that the pure transport equation (9), the one which does posses the frame-invariant property, can be viewed as a limiting equation by letting γ → 0 in (16). The latter could be viewed as a (PDE based) regularization of the pure transport system, which is inconsistent to the classical description, such as those with sharp interface formulations. In the case of equal density, the last term on the right-hand side vanishes and hence we recover the usual ACNS model (e.g., [25]).
We consider the usual no-slip boundary condition for velocity and natural boundary condition for phase variable Summing up, the convective Allen-Cahn equation (16) and the momentum equation (15) with the boundary condition (17), together with the incompressibility constraint form a complete system for (u, p, φ) with given function ρ(φ) satisfying certain physically consistent conditions (e.g. (24)-(25)).

Remark 2.
In calculations of (13), we take integration by parts first and then change the coordinates. Alternatively, we can change the coordinates first and then integrate by parts which leads to the following: Then it results in the following macroscopic force balance equation: wherep is due to the Lagrangian multiplier. More interestingly, if we use both ways to one half of the first term on the second line of (13), we obtain the following momentum equation: Then the system consisting of (20), (18) and the usual convective Allen-Cahn equa- is also thermodynamically consistent. Indeed, this system enjoys the same energy law as (1)-(3). In a future work, we will study this system under a mass-preserving constraint, i.e., (21) with presence of a Lagrangian multiplier ξ(t) due to a nonlocal constraint Ω φdx = Ω φ 0 dx.

Existence of weak solutions.
We now aim to investigate theoretically the ACNS system derived in the previous section. In this section, we prove existence of a global weak solution to problem (1)-(6).
3.1. Functional setting and basic assumptions. First, we give some notions and notation. In the sequel, Ω denotes a bounded domain with smooth boundary in R n (n = 2, 3). L p (Ω) and W m,p (Ω) with m being an integer and p > 1 denote the usual Lebesgue spaces and Sobolev spaces, respectively and as usual, W m,2 will be denoted by H m . With no ambiguity, the corresponding vector Lebesgue and Sobolev spaces will also be denoted by L p and W m,p . On the other hand, we consider the strictly positive operator A : a.e. on ∂Ω} and for any s ∈ R, D(A s ) is well defined (see [40]). Moreover, we will use the fact that

Now, we introduce the Hilbert spaces
with scalar product and the norm of L 2 (Ω) and endowed with the product and norm Let the operator P be the Leray-Helmholtz projection of L 2 (Ω) onto H and consider the eigenvalue problem: It is well-known that −P ∆ is a positive, self-adjoint operator in H and there exist two sequences {λ j } j=1,2... and {ω j } j=1,2... such that, for every j ≥ 1, λ j ≥ 0 is an eigenvalue and ω j is a corresponding eigenfunction, λ j is nondecreasing, tending to infinity as j → +∞ and {ω j } is orthonormal and complete in H. Now we introduce the assumptions on the parameter functions. Without loss of generality, we choose the typical double-well potential for the free energy as follows: with some ε > 0 and for the average density function, we assume that with ρ 1 < ρ 2 being two positive constants and the following exterior convexity: At last, all the physical parameters in the system are normed to be one for the sake of simplicity.
Before study the existence of weak solutions to problem (1)-(6), we first state an auxiliary result concerning with a convective Allen-Cahn equation. For the reader's convenience, its proof is reported in the Appendix.
Moreover, it holds that Now let us formulate the result concerning with the solvability of problem (1)-(6) in the sense of weak solutions.
in Ω, problem (1)-(6) admits at least one solution in the sense of Definition 3.1.
The proof of Theorem 3.3 consists of several steps. First, we construct an approximation problem via Galerkin scheme and we obtain an approximation solution by a fixed point argument. Then, we prove the compactness of the approximation solutions. Last, we obtain a weak solution by passing to the limit.
Step 1. Let m > 0 be a fixed integer and set u m ( Proof. We use a fixed point argument. Denote with T 0 , M > 0 to be specified later. Suppose that we are given h m : We then look for the solution (u m , φ m ) that satisfies the following (nonlinear) auxiliary problem: First, noticing that ω i is regular since ∂Ω is smooth, we have v m ∈ C([0, T 0 ]; C 2 (Ω)). Then, by Lemma 3.2, we get a solution φ m (t) ∈ C([0, T 0 ]; H 1 (Ω))∩L 2 (0, T 0 ; H 2 (Ω)) satisfying the second equation of (33) and moreover, we have |φ m (x, t)| ≤ 1 pointwise in Ω × [0, T 0 ]. As a result, since ρ ∈ C 1 (R), there exists a positive constant β Moreover, multiplying the second equation in (33) by ∂ t φ m and integrating by parts, one obtains d dt Ω where C 1 is a constant depending on max 1≤i≤m ω i L ∞ (Ω) and β. Then it follows from the Gronwall inequality that Concerning the solvability of (37), we may first state the following result which can be proved in the same way as in [22]. (38) is nonsingular and each component of its inverse belongs to C[0, T 0 ].

Then noticing that
Moreover, we have g m (t) ∈ U m with proper M . Indeed, multiplying (37) by g i m (t)(A −1 ) ij and taking summation in i and j yields that Since (A −1 ) ij , B j kl are bounded from both above and below, we obtain by Young's inequality that where C 2 depends on m, β and ω i C 1 (Ω) . Then, applying the Gronwall inequality again and using (36), we get for any t ∈ [0, In this way, we can define a nonlinear map L from U m (T m ) into itself such that L(h m ) = g m . Moreover, it is not difficult to see that L maps U m (T m ) into a bounded set in (W 1,p (0, T m )) m for any p > 1. Therefore, L is compact due to the Sobolev embedding theorem. Then by the fixed point theory of Leray and Schauder, we infer that there is a fixed point of L in U m (T m ) which means problem (AP) has a local solution. This completes the proof of Proposition 1.
Next, we prove that the local solution of problem (AP) can be actually extended to [0, T ] for any T > 0. To this aim, we need to get some a priori estimates. First, in the same way as before, Then, multiplying the first equation in (32) by g i m and second equation in (32) by ∂ t φ m + u m · ∇φ m , respectively, then integrating by parts and adding the resultant up, we get It follows by an integration with respect to time that Observing that and ρ(φ m ) is bounded from above and below, we infer that where C depends only on φ 0 , u 0 and Ω. Hence, the local solution of problem (AP) can be extended globally.
Step 2. With the uniform estimate (41) at hand, observing the 3D Ladyzhenskaya inequality and the Gagliardo-Nirenberg inequality and for any p > 1, Therefore, we infer that there is a subsequence of (u m , φ m ), still denoted by (u m , φ m ), such that and weakly in L Moreover, and due to the Aubin-Lions lemma, we infer that and hence φ m , ∇φ m → φ, ∇φ, a.e. in (0, T ) × Ω.
Next, we prove the following lemma concerning uniform estimates of translations (cf. [9]) which plays an important role in deriving the strong convergence of u m . Lemma 3.5. For arbitrary δ ∈ (0, T ), it holds that where C is independent of m and δ.
Proof. Denote ρ m = ρ(φ m ) for the sake of simplicity. For fixed δ and t such that 0 < δ < T and 0 ≤ t ≤ T − δ, we have we get Observing that and making use of the property 1 β ≤ ρ m ≤ β, we derive that Integrating (48) with respect to t from 0 to T − δ, we claim that for any J i on the right-hand side of the above inequality, the following estimates hold: where constants κ i are independent of m and δ. First, by Hölder's inequality, we have Hence, we have .
On the other hand, noticing that by Hölder's inequality, In an analogous way, observing that it follows by change the order of integration that Finally, for J 4 , noticing that This completes the proof. Now we may apply the following compactness result which was developed by J.L. Lions [24] to establish existence results for incompressible fluids with non-constant density.
Thanks to Lemma 3.5 and Lemma 3.6, we conclude that u m is compact in L 3 (0, T ; L 9 4 (Ω)). Therefore, we can extract a subsequence, still denoted by u m , such that and hence u m → u a.e. in (0, T ) × Ω.
Step 3. In order to examine the weak limits of the nonlinear terms in (AP), we will make use of the following lemma.
Lemma 3.7. Suppose that Ω is a bounded domain in R n . Let {w k (x)} ∞ k=1 be a bounded sequence in L p (Ω) (1 ≤ p < ∞) such that w k almost everywhere converges to w. Then w also belongs L p (Ω) and w k weakly converges in L p (Ω) to w. If p = ∞, the conclusion becomes that w k weak star converges to w.
Applying this lemma and in view of (43)-(46) and (50), we may easily deduce that
Note that holds for any fixed j, ω j ∈ C 2 (Ω). For the first term on the right-hand side of (52), we deduce by (50) that For the second term of (52), due to Hölder's inequality, we have where by Lebesgue's dominated convergence theorem, As a result, On the other hand, since ρ (φ)uηω j ∈ L 3 (0, T ; L 9 4 (Ω)), it follows from (44) that Now, based on the above argument, one easily get by passing to the limit in (51) that holds for v = ω 1 , ω 2 , .... Since {ω j } ∞ j=1 is the basis of V and each term in (54) is valid for v ∈ V, we conclude that (54) is still true for any v ∈ V. On the other hand, passing to the limit in the second equation of (32), we easily conclude that (29) holds in L   The existence of a strong solution can be proven in a similar procedure via compactness argument as done in the previous section. Hence, we only need to show some necessary higher order a priori estimates and limiting passage.
First, multiplying the first equation of (32) by (d/dt)g j m (t) and λ j g j m (t), > 0 and taking summation over j = 1, ...m, we obtain that We observe that By the 2-D Galiardo-Nirenberg inequality u L ∞ (Ω) ≤ ∆u It follows from (42) and (56) that with δ > 0 to be chosen later. On the other hand, we notice that where D i denotes the i−th order derivatives with respect to x and similarly, For the last term on the right-hand side of above two inequalities, we have where we use the 2D Gagliardo-Nirenberg inequalities: and global boundedness (41). Now, multiplying the second equation in (32) by −∂ t ∆φ m , integrating by parts and applying Young's inequality, we obtain that For the last term of the above inequality, we observe that and recalling (57), Noticing that (58) and on the other hand, from the second equation of (32) and picking δ ∈ (0, 1 4C0 ) in above estimates, we infer that Collecting all above inequalities, picking ∈ (0, 1 4β ) and δ sufficiently small such that δ ≤ 8C0(β 2 +β+1) , we finally arrive at we deduce from the Gronwall inequality that Therefore, recalling (41), we may deduce that ∂ t u m ∈ L 2 (0, T ; L 2 (Ω)) and are uniformly bounded. It follows that φ m → φ weakly* in L ∞ (0, T ; H 2 ) and weakly in L 2 (0, T ; H 3 ) ∂ t φ m → φ t weakly* in L ∞ (0, T ; L 2 ) and weakly in L 2 (0, T ; H 1 ) u m → u weakly* in L ∞ (0, T ; V) and weakly in L 2 (0, T ; H 2 ) ∂ t u m → u t weakly in L 2 (0, T ; L 2 (Ω)).
(61) Moreover, by Aubin-Lions lemma we may infer that and u m → u strongly in L 2 (0, T ; V). Hence, it is not difficult to obtain that By Lebesuge convergence theorem, we have as m → +∞, On the other hand, since ρ(φ)ω j η(t) ∈ L 2 (0, T ; L 2 (Ω)), it follows from (61) that Finally, passing to the limit in (32), based on above arguments, we may conclude that (u, φ) is a strong solution in the sense of Definition 4.1. This completes the proof of Theorem 4.2.
In addition, under the assumption that ρ is local Lipchitz, we may prove the continuous dependence of solutions on the initial data.
in Ω. In addition, assume that there exists a generic constantC such that Then for t ∈ [0, T ], it holds that where the constant L may depend on u i0 V , φ i0 H 2 , Ω,C.
In a similarly way as before, we obtain that We observe that H 2 ), H 2 ∇u 2 and In view of assumption (62) and (56), we deduce that and On the other hand, multiplying (63) by φ and integrating by parts yields that Therefore, denoting we may infer from all above inequalities that where with C being a positive constant depending on u i0 , φ i0 . Then, it follows from Gronwall's inequality that y(t) ≤ e Lt y(0).
Consequently, we conclude that This completes the proof. Remark 4. It is our future interest to study whether system (1)-(3) provides an applicable approximation to obtain existence of solutions to (1)-(2) coupling with the pure transport equation (9), since (3) can be viewed as a relaxation of (9).

5.
Longtime behavior of the strong solutions. The study of asymptotic behavior of the solution to nonlinear evolution equation as time goes to infinity can be roughly divided into two categories. The first category is to investigate the asymptotic behavior of all solutions when initial data vary in any bounded set. The second category is to investigate the asymptotic behavior of single trajectory starting from any point in the space. In this section, we will consider the longtime behavior of the global strong solutions in the 2D case from both points of view.

Existence of maximal attractors. Let
which is a complete metric space when equipped with the metric induced from the norm of V × H 2 . Then by Theorem 4.2 and Theorem 4.3, we conclude that the nonlinear semigroup S(t) defined by problem (1)-(6) maps X into itself and is strongly continuous. Moreover, we show that S(t) possesses an absorbing set in X.

Proposition 2.
There is a bounded set B 0 in X such that B 0 is absorbing in X, i.e., for any bounded set B in X, there exists t 1 > 0 depending on B such that On the other hand, multiplying (3) by φ and integrating over Ω, we obtain that From here and below, K i will denote generic constants that depend on Ω, β and the coefficients of the system, but are independent of the initial data. Now, letting with 1 > 0 to be chosen later. Then, we deduce from above that By Poincaré's inequality, ∇u 2 ≥ λ 1 u 2 , we infer that On the other hand, noticing that and picking 1 = min{ λ1 4K1 , 1 2 }, we arrive at with δ 1 = min{ λ1 β , 2 1 } and K 2 = 4 1 |Ω|. Hence, we get Λ 1 (t) ≤ Λ 1 (0)e −δ1t + K 2 . Since we infer that there exists t 0 > 0 depending on Λ 1 (0) such that when t ≥ t 0 , Moreover, noticing that In the same way as deriving (60), we may deduce that for t ≥ t 0 Then making use of the uniform Gronwall inequality, we get which together with (71) implies that when t ≥ t 1 := t 0 + 1 The proof is finished.
Proposition 3. The operator S(t) is uniformly compact for t large, i.e., for every bounded set B contained in X, there is t 2 > 0 that may depend on B such that t≥t2 S(t)B is relatively compact in X. Proof. Since we are only interested in the behavior of solutions after a sufficiently large time, without affecting the analysis, we may assume that (72) holds for t ≥ t 1 = 0. Now multiplying (1) by −∆u t and integrating by parts, keeping the uniform estimates (72) in mind, we obtain that Here and in the sequel, C denotes a constant that may vary from line to line and depend on the initial data but is independent of T . We observe by the 2D Agmon inequality, (58) and (72) that and by 2D Gargliardo-Nirenberg inequality and Hölder's inequality that and Ω |u| 2 |D 2 u| 2 dx ≤ C ∆u 4 .
Recalling that (see (59)) Now, keeping (72) in mind and applying the uniform Gronwall inequality again, we infer that for t ≥ t 2 := t 1 + 1, it holds that and hence from the equations we have, for t ≥ t 2 := t 1 + 1,  (1)- (6) possesses in X a maximal attractor A which is compact.

Convergence to equilibria.
In this part, we analyze the asymptotic behavior of single trajectory as time goes to infinity. First, we define the ω−limit set as follows.
Next, we introduce the following stationary Allen-Cahn problem −∆ψ + f (ψ) = 0 ∂ ν ψ| ∂Ω = 0 (78) and for any ψ ∈ H 2 N (Ω), we define its associated functional as Then, it is not difficult to obtain the following lemma by calculus of variations and standard regularity theory for elliptic equations.
Lemma 5.2. Suppose that ψ ∈ H 2 (Ω) satisfies (78). Then ψ is a critical point of the functional E over H 2 (Ω). Conversely, if ψ ∈ H 2 (Ω) is a critical point of E and ∂Ω is of class C ∞ , then ψ ∈ C ∞ (Ω) and it is a classical solution of problem (78).
Remark 5. The decay property of the velocity u as time tends to infinity can be obtained by the energy method (see the proof below). However, convergence for the phase field parameter φ is usually nontrivial because the structure of the set of equilibria may be complicated and the solutions to elliptic problem like (78) may form a continuum if the spatial dimension n ≥ 2 (see e.g., [21,Remark 2.3.13]). Since our problem enjoys a dissipative energy identity, we can achieve the goal via the so-called Lojasiewicz-Simon approach (cf. e.g., [18,19,38]). One advantage of this approach is that we can obtain the convergence result without investigating the structure of equilibria.
The version of the Lojasiewicz-Simon inequality we need here is given by the following lemma whose proof can be found in [40] in detail. Hence we omit its proof here. We remark that in general, the lemma is valid when the nonlinear term f (·) is real analytic that is obviously satisfied by our specific potential (23).
Lemma 5.5. Let ψ ∈ H 2 N (Ω) be a critical point of E. Then there exists constants θ ∈ (0, 1 2 ) and σ > 0 depending on ψ such that, for any φ ∈ H 2 N (Ω), when Proof of Theorem 5.4. First, we show that u(t) → 0 in V as time goes to infinity via energy method. Observe that by a direct computation and recall (76)-(77), d dt On the other hand, integrating the energy identity (69) with respect to time, we get ∞ 0 ∇u 2 dτ ≤ C.
(84) Since for any strong solution (u(t), φ(t)), total energy E(u(t), φ(t)) is non-increasing in t and E is bounded below, we may infer that E(u(t), φ(t)) → E ∞ for some constant E ∞ which should equal the constant value of E on ω(u 0 , φ 0 ) because of the definition of ω−limit set and the continuity of E in X. Thus, for any (0, ψ) ∈ ω(u 0 , φ 0 ), it holds that E(ψ) = E(0, ψ) = E ∞ (85) and moreover, from (84) and uniform upper bound of ρ, we easily deduce that Due to Lemma 5.3-5.5, for any (0, ψ) ∈ ω(u 0 , φ 0 ), there exists a σ ψ > 0 and θ ψ ∈ (0, 1 2 ) such that the inequality (83) holds for The union of balls B 1 (0) × B σ ψ (ψ) forms an open cover of ω(u 0 , φ 0 ), where B 1 (0) denotes the unit ball centered at origin in V. Due to the compactness of ω(u 0 , φ 0 ) in V × H 2 (Ω), we can find a finite sub-cover where the constants σ i , θ i corresponding to ψ i in Lemma 5.5 are indexed by i. From the definition of ω(u 0 , φ 0 ), there is a sufficiently large t 0 > t 2 such that Now, taking θ = min m i=1 θ i ∈ (0, 1 2 ), using Lemma 5.5 and (85), we deduce that for all t ≥ t 0 , Observe that d dt E(u, φ) + ∇u 2 + φ t + u · ∇φ 2 = 0 (89) and from the equation for φ, where we use the 2D Ladyzhenskaya inequality and the uniform bounds for u . We infer that there exists some constant ζ > 0 such that Let Then it follows from (91) and (86) that for any t ≥ t 0 , On the other hand, recalling (88), we deduce that for t ≥ t 0 where we make use of 1 < 1 1−θ < 2, the uniform pointwise bounds of ρ(·) and the uniformly bounds of u in H together with Poincaré's inequality. It follows from (92)-(93) that ∞ t0 Z 2 (s)ds ≤ CZ Now applying [17,Lemma 7.1], we conclude that Z(t) ∈ L 1 (t 0 , +∞). As a result, using the 2D Ladyzhenskaya inequality and Poincaré's inequality again, we obtain that It follows that for any ψ such that (0, ψ) ∈ ω(u 0 , φ 0 ), it holds that This means that the whole trajectory converges to ψ in L 2 . Due to the uniqueness of the limit, we conclude that It remains to prove the estimate on convergence rate. To this aim, we first notice that due to (91) and (93), it holds that which implies that Furthermore, a direct computation together with (93) yields Then an integration with respect to time of the above inequality together with (94) gives Adjusting the constant C properly, we may obtain that Now, letting ϕ = φ − ψ, we have Note that for any v ∈ V, Then multiplying (97) by ϕ t + u · ∇φ and taking the test function of (55) as u, then adding the resultants together, we obtain that +C Ω |ϕ| 2 |φ 2 + φψ + ψ 2 | 2 + |∇∆ψ| 2 + 1 Observing that from (97), uniform boundedness of u in H and pointwise boundedness for φ and ψ, it holds Then from (99) together with (98), we arrive at where C depends on ψ C 3 . Now, letting Y(t) = Ω ρ(φ)|u| 2 + |∇ϕ| 2 and noticing that Y(t) ≤ β u 2 + ∇ϕ 2 ≤ C( ∇u 2 + ϕ 2 H 2 ), we derive that there exists some λ 0 > 0 such that Adjusting the constant C properly, we get for t ≥ 0 Consequently, Once obtaining the above convergence rate of lower order, we are allowed to further derive the same convergence rate for higher order norms using energy method in an analogous way as in previous sections and hence we omit the details here (cf. [19,38] Appendix. For the reader's convenience, we report a proof for Lemma 3.2. Proof. For any 0 < < 1 2 , let J be the standard mollifier and ρ = J * ρ. Then it is easy to check that ρ (s)s ≥ 0, |s| ≥ 1 + (100) and ρ → ρ in C 1 ([−2, 2]).
Moreover, we may construct a smooth sequence φ 0 ∈ H 2 N (Ω) such that lim with δ > 0 to be chosen later. For any ψ ∈ W δ , consider the following auxiliary linear parabolic equation: Since ρ ∈ C 3 (R), it is easy to check that g ∈ C([0, δ]; D(A)). Denoting {e −tA } t≥0 the semigroup generated by A, we have the following expression of φ : φ (t) = e −tA φ 0 +