Finite Volume approximation of a two-phase two fluxes degenerate Cahn-Hilliard model

We study a time implicit Finite Volume scheme for degenerate Cahn-Hilliard model proposed in [W. E and P. Palffy-Muhoray. Phys. Rev. E, 55:R3844-R3846, 1997] and studied mathematically by the authors in [C. Canc\`es, D. Matthes, and F. Nabet. Arch. Ration. Mech. Anal., 233(2):837-866, 2019]. The scheme is shown to preserve the key properties of the continuous model, namely mass conservation, positivity of the concentrations, the decay of the energy and the control of the entropy dissipation rate. This allows to establish the existence of a solution to the nonlinear algebraic system corresponding to the scheme. Further, we show thanks to compactness arguments that the approximate solution converges towards a weak solution of the continuous problems as the discretization parameters tend to $0$. Numerical results illustrate the behavior of the numerical model.


The continuous model
We consider a mixture made of two incompressible phases evolving in a bounded and connected polygonal open subset Ω of R 2 and on a time interval [0, ], where is an arbitrary finite time horizon. The composition of the fluid is described by the volume fractions = ( 1 , 2 ) of the two phases. Since the whole volume Ω is occupied by the two phases, the following constraint on the holds In the above equation, > 0 denotes the viscosity of the phase , is its chemical potential (which is one of the unknown of the problem), while Ψ ∈ 1 (Ω) is a given external potential acting on phase that is assumed be independent on time for simplicity. For Ψ , one can typically think about gravity, that is Ψ ( ) = − g · with the density of phase and g the gravitational vector. The coefficient ≥ 0 is a given parameter quantifying the thermal agitation of phase . The limit case = 0 is called the deep-quench limit in the Cahn-Hilliard literature.
The difference of the phase chemical potentials is given by the following expression where > 0 and > 0 are given coefficients governing the characteristic size of the transition layers between patches of pure phases { 1 = 0} and { 1 = 1}. Typically, is assumed to be small in comparison to . Equation Finally to close the system, we impose an initial condition 0 = ( 0 1 , 0 2 ) on the volume fractions by setting The initial profiles 0 ∈ 1 (Ω) are assumed to be nonnegative with 0 1 + 0 2 = 1 in Ω, and we assume that both phases are present at initial time, i.e.,

Fundamental estimates and weak solutions
As a preliminary to the study of the numerical scheme, we derive formally at the continuous level some a priori estimates. Their transposition at the discrete level will be key in the numerical analysis to be proposed in what follows. Equation (1.2) can be rewritten under the form + ∇ · = 0, with = − ∇ ( + Ψ + log( )) . ∫︁ Ω |∇ ( + Ψ + log( ))| 2 d ≥ 0.
The second term in the right-hand side can be rewritten as while we can make use of (1.3) to rewrite the first term as Therefore, we obtain the energy/energy dissipation relation d d E( ) + D( , ) = 0 for all ≥ 0, (1.11) where the energy functional E( ) is defined by (1.12) A straightforward consequence of (1.11) is that ↦ → E( ( )) is non-increasing along time, and thus that We deduce from previous inequality that the energy is bounded, hence a ∞ ((0, ); 1 (Ω)) estimate on . The energy/energy dissipation estimate (1.11) is not sufficient to carry out our mathematical study since it only provides a weighted estimate on the chemical potentials ∑︁ ∈{1,2}
The first term is in 2 ((0, ); 1 (Ω)) as the product of an element of 2 ( ) with an element of ∞ ((0, ); 2 (Ω)), while the second term is in 2 ( ) since 0 ≤ ≤ 1 and thanks to (1.14). As a consequence, ∇ is bounded in 2 ((0, ); 1 (Ω)). Making use of the Poincaré-Sobolev estimate (recall that has zero mean for all time, cf. (1.6), and that Ω ⊂ R 2 ), we obtain that is bounded in 2 ( ). To get the desired 2 ( ) estimate on 1 , it only remains to check that belongs to 2 ( ) thanks to the 2 ( ) estimates on and 1 − 2 together with 0 ≤ 2 ≤ 1. The interest of the above formal calculations is twofold. First, our scheme has been designed so that all these calculations can be transposed to the discrete setting. The corresponding a priori estimates will be at the basis of the numerical analysis proposed in this paper. Second, these estimates provide enough regularity on the solution to give a proper notion of weak solution to the problem.
as well as ∫︁ ∫︁ The existence of a weak solution has been established in [8] by showing the convergence of a minimizing movement schemeà la Jordan et al. [23]. Note that in [8], the case of a convex three-dimensional domain Ω is also addressed, but it relies on the fact that the 2 ( ) estimate on ∆ 1 yields a 2 ((0, ); 2 (Ω)) estimate on for which we don't have an equivalent at the discrete level. This is why we restrict our attention on the case Ω ⊂ R 2 (but not necessarily convex) in this paper.

Some words about the model
Before entering the core of the paper, which is devoted to the convergence analysis of a finite volume scheme, let us briefly discuss the model under consideration, and in particular its difference with respect to the usual Cahn-Hilliard system. We refer to [8,29] for a more developed discussions on this purpose. The classical degenerate Cahn-Hilliard equation which is the closest one to our system (1.1)-(1.7) writes (1.18) where the degenerate mobility is given by ( ) = , and the generalized chemical potential This system has to be completed with one initial condition on 0 1 and boundary conditions The existence of a solution to this problem has been addressed in [14,25]. Let us come back to our system (1.1)-(1.7). Denote the total flux by tot = 1 + 2 , then summing (1.9) over ∈ {1, 2} yields ∇ · tot = − ( 1 + 2 ) = 0 in (0, ) × Ω (1.20) owing to (1.1). After some elementary calculations, the conservation (1.2) for the phase 1 rewrites where ( ) = 2 1+( 2 − 1) . The equation (1.21) differs from (1.18) by the addition of a nonlinear transport term driven by a divergence free vector field. Both systems can be reinterpreted as Wasserstein-type gradient flows [1] of the energy E( ) for different geometries: -the Wasserstein distance with quadratic cost with the constraint tot = 0 for the classical Cahn-Hilliard system, cf. [25]; -the Wasserstein distance with quadratic cost with the less stringent constraint ∇ · tot = 0 for the system (1.1)-(1.7), cf. [8,29].
The additional degree of freedom tot allows the energy E( ) to decrease faster along the trajectories, as highlighted in [8]. Finally, let us point the recent contribution [7] where the convergence of a minimizing movement scheme is addressed for a closely related model where the Cahn-Hilliard energy is replaced by the singular de Gennes-Flory-Higgins energy.

Finite volume approximation and main results
Prior to presenting the scheme and stating our main results, that are the existence of a discrete solution to the scheme and the convergence of the corresponding approximate solutions towards a weak solution to the problem (1.1)-(1.7), we introduce some notations and requirements concerning the mesh.

(Super)-admissible mesh of Ω and time discretization
Let us first give a definition of what we call an admissible mesh.
such that the following conditions are fulfilled.
(i) Each control volume (or cell) ∈ is non-empty, open, polygonal and convex. We assume that We denote by the 2-dimensional Lebesgue measure of . (ii) Each edge ∈ ℰ is closed and is contained in a hyperplane of R 2 , with positive 1-dimensional Hausdorff (or Lebesgue) measure denoted by = ℋ 1 ( ) > 0. We assume that ℋ 1 ( ∩ ′ ) = 0 for , ′ ∈ ℰ unless ′ = . For all ∈ , we assume that there exists a subset ℰ of ℰ such that = ⋃︀ ∈ℰ . Moreover, we suppose that ⋃︀ ∈ ℰ = ℰ. Given two distinct control volumes , ∈ , the intersection ∩ either reduces to a single edge ∈ ℰ denoted by | , or its 1-dimensional Hausdorff measure is 0.
(iii) The cell-centers ( ) ∈ are pairwise distinct with ∈ , and are such that, if , ∈ share an edge | , then the vector = − | − | is orthogonal to | . (iv) Given two cells , ∈ sharing an edge = | , we assume that the straight line joining and crosses the edge in its midpoint .
Let us introduce some additional notations, some of them being depicted on Figure 1. The size of the mesh (which is intended to tend to 0 in the convergence proof) is defined by ℎ = max ∈ ℎ , with ℎ = diam( ). Given two neighboring cells , ∈ sharing an edge = | , we denote by = | − | whereas = | − | ≤ . The transmissivities and of the edge are respectively defined by = and = . The diamond and half diamond cells are defined as the convex hulls of { , , } and { , } respectively. Denoting by (resp. ) the 2-dimensional Lebesgue measure of (resp. ), we will use many time the following elementary geometric properties: = 2 and = 2 . We also denote by ℰ ,int the subset of ℰ made of the internal edges such that there exists ∈ such that = | , and by ℰ int = ⋃︀ ∈ ℰ ,int . Even though this is absolutely not necessary, we choose to restrict our attention to the case of uniform time discretizations in the mathematical proofs in order to reduce the amount of notations. In what follows, we set ∆ = / and = ∆ for ∈ {0, . . . , }. The integer is intended to be large and even to tend to +∞ in the convergence proof.
Remark 2.2. Condition (iv) above enforces an additional restriction with respect to the classical definition of finite volumes meshes with orthogonality condition (iii). Meshes satisfying this condition in addition to the more classical assumptions (i)-(iii) is called super-admissible following the terminology introduced in [17]. It is for instance satisfied by cartesian grids or by acute triangulations. However, (iv) is in general not satisfied by Voronoï meshes. This condition appears for technical reasons related to the construction of a strongly convergent SUSHI discrete gradient [17], see Proposition 4.2 later on. On the other hand, this condition was recently pushed forward in [20] to show the consistency of the discrete optimal transportation geometry [26] hidden behind our work with the continuous optimal transportation geometry [31] in which our system (1.1)-(1.7) has a gradient flow structure [8].

A two-point flux approximation finite volume scheme
The scheme we propose is a cell-centered scheme based on two-point flux approximation (TPFA) finite volumes. At each time step ∈ {1, . . . , }, then unknowns are located at the centers of the cells ∈ .
Given discrete volume fractions −1 = (︁  As highlighted in the formal calculations presented in Section 1.2, the analysis requires the use of the logarithm of the volume fractions. To this end, the volume fractions , have to be strictly positive for ≥ 1. To ensure this property, some thermal diffusion is needed, see Lemma 3.2. In the case where = 0, then one needs to introduce a small amount of numerical diffusion by setting where > 0 is a parameter that can be fixed by the user. Equation (1.2) is then discretized into for all ∈ and ∈ {1, 2}. In the above relation, we used the following discretization of the external potential: Edge values , of the discrete volume fractions also appear in (2.3). Rather than using upstream values of the volume fractions as in our previous work [9], we make use of a logarithmic mean, i.e., for all = | ∈ ℰ int and ∈ {1, 2}. This particular choice of the edge value is closely related to the one introduced in the early work [22] for the approximation of the thin film equation, and fits with the one suggested in [26,28] and used in a closely related context to ours in [27].
Note that the repulsive term (second in the right-hand side) is discretized in an explicit way for stability issues that will appear clearly later on. The constraint (1.1) is discretized in a straightforward way by imposing The last equation to be transposed in the discrete setting is (1.6), which is translated into

Main results
Before addressing the convergence of the scheme, we focus first on the case of a fixed mesh and time discretization. The scheme (2.3)-(2.7) yields a nonlinear system on ( , ). The existence of a solution to this nonlinear system is far from being obvious. The existence of such a solution and some important properties of the discrete solution mimicking the properties highlighted in Section 1.2 are gathered in the first theorem of this paper. (ii) positivity: (iii) energy decay: where the discrete energy E ( ) is defined by where is given by (1.10), and the discrete dissipation is defined by The existence of a solution to the scheme for each time step allows to reconstruct an approximate solution ( ,Δ , ,Δ ) with ,Δ = ( 1, ,Δ , 2, ,Δ ) and ,Δ = ( 1, ,Δ , 2, ,Δ ) by setting The approximate solutions are expected to approximate the continuous solution to (1.1)-(1.7). Our second theorem gives a mathematical foundation to this statement. It requires the introduction of a suitable sequence of discretizations of . In what follows, we denote by (︀ , ℰ , ( ) ∈ )︀ ≥1 a sequence of admissible meshes of Ω is the sense of Definition 2.1. We assume that ℎ tends to 0 as → ∞ as well as the following regularity requirements: -shape regularity of the cells: there exists a finite > 1 such that and such that -boundedness of the number of edges per element: there exists ℓ ⋆ ≥ 3 such that 14) The combination of a sequence (︀ , ℰ , ( ) ∈ )︀ ≥1 fulfilling (2.11)-(2.14) together with a time step ∆ = / and → +∞ as → ∞ is said to be a regular discretization of if it moreover satisfies the inverse CFL condition (3.18). The convergence properties stated in Theorem 2.4 are weaker than what is practically proved in the paper. The statement of optimal convergence properties would require the introduction of additional material that we postpone to the proof in order to optimize the readability of the paper. The remaining of the paper is organized as follows. In Section 3, we work at fixed mesh and time step. We derive some a priori estimates and show the existence of (at least) one solution to the scheme thanks to a topological degree argument. Next in Section 4, we show thanks to compactness arguments that the approximate solution converge towards a weak solution to the scheme. Finally, we present in Section 5 some numerical simulations.

A PRIORI estimates and existence of a discrete solution
In Section 3.1, we first derive some a priori estimates on the solutions to the scheme (2.3)-(2.7). These estimates will be at the basis of the existence proof for a discrete solution to the scheme in Section 3.2, but also of the convergence proof carried out in Section 4.

A priori estimates
This section is devoted to the derivation to all the a priori estimates needed in the numerical analysis of the scheme. The first of them is the global conservation of mass, which is a consequence of the local conservativity of the scheme.
Proof. Summing (2.3) over ∈ and using the conservativity of the scheme leads to A straightforward induction and the definition (2.1) of 0 , then provides (3.1).
Our second lemma shows that the volume fractions are positive.  .7), then > 0 (the later having been established for all ≥ 1 in Lem. 3.1), and suppose for contradiction that , = min ∈ , ≤ 0, so that , = 0 for all ∈ ℰ ,int . Therefore, on this specific control volume , the scheme (2.3) reduces to The left-hand side is nonpositive, and even negative unless , = , ≤ 0 for all the neighbouring cells of . We can thus iterate the argument and show that , ≤0 for all ∈ , which provides a contradiction with the property ∑︀ , > 0 established in Lemma 3.1.
As a consequence of Lemma 3.2, the quantities log (︀ , )︀ have a sense. They will be used many times along the paper. Our next lemma consists in discrete counterparts of the energy/energy dissipation relations (1.11)-(1.13).
where the discrete energy E and the discrete dissipation D are defined by (2.8) and (2.9) respectively.
The boundedness of the discrete energy E ( ) provides a discrete ∞ ((0, ); 1 (Ω)) estimate on the volume fractions, as established in the next corollary.
Proof. As a straightforward consequence of Lemma 3.3, the energy is decaying along the time steps, so that Owing to Lemma 9.4 of [16], there exists 2 depending only on such that whereas Jensen's inequality ensures that Let us now focus on the quantification of the production of mixing entropy at the discrete level. Our next lemma provides a discrete counterpart to Estimate (1.15).
Lemma 3.5. There exists 3 depending only on Ω, , , , Ψ , , , and 0 1 such that As a consequence, there exists 4 depending only on Ω, , , , Ψ , , , and 0 1 such that where we have set It follows from the convexity of that . (3.14) The particular choice (2.4) for , was fixed so that Therefore, using (2.6) and Cauchy-Schwarz inequality, we deduce that where the last inequality is a consequence of Corollary 3.4 and of estimate which itself is a consequence of Lemma 9.4 from [16] and of the 1 (Ω) regularity of the external potentials Ψ . Similarly, one can rewrite Thanks to the relation (2.5), it turns to Using the fact that 0 ≤ −1 1, The combination of (3.13)-(3.17) in (3.12) provides (3.10). Let us now focus on estimate (3.11). Equality (2.5) gives Since 0 ≤ −1 1, ≤ 1 and the logarithmic function is increasing, estimate (3.10) concludes the proof.
The following lemma is a transposition to the discrete setting of the weighted estimate (1.14) on the chemical potentials.
Lemma 3.6. There exists 5 depending only on , , 0 , Ψ , , Ω, and such that Proof. Definition (2.9) of D ( , ) together with inequality ( + + ) 2 ≤ 3( 2 + 2 + 2 ) yield Owing to Lemma 3.3, the first term of the right-hand side is bounded by which is bounded as already seen in the proof of Corollary 3.4. On the other hand, since 0 ≤ , ≤ 1, one has Finally, thanks to Lemma 3.5.
Relation (2.6) guarantees that the sum of the volume fractions is constant equal to 1 in the cells. But this is no longer true on the edges. As shown in the following lemma, the sum of the edge volume fractions is always lower or equal to 1. Assume for instance that for some = | , 1, = 1 and 1, = 0, then both 1, and 2, are equal to 0. This degeneracy may lead to severe difficulties in the effective resolution of the nonlinear system provided by the scheme. Next lemma shows that this situation can not be encountered provided the time step is large enough with respect to the size of the mesh. The estimate we provide is based on the worst case scenario and is thus extremely pessimistic. Practically, the inverse CFL condition (3.18) is not needed as soon as the ratio / is large enough with respect to the size of the discretization. Lemma 3.7. Assume that there exists > 1 such that then there exists ∈ (0, 1) depending on ⋆ , ⋆ , ℓ ⋆ and such that As a consequence, there exists ⋆ > 0 depending only on such that

20)
Proof. Let us first establish (3.19). As a consequence of Lemma 3.5, there holds Plugging it in (3.21) and using (2.14) yields Let us now turn to the proof of (3.20). If , = , = , , we have immediately 1, + 2, = 1. Otherwise, the inequality 1, + 2, ≤ 1 follows directly from the fact that the logarithmic mean is smaller than the arithmetic one. Define the continuous function : Thus it remains bounded away from 0 by some ⋆ depending only on . Then (3.20) follows from (3.19) and (3.22).
With Lemma 3.7 at hand, we are in position to prove our next lemma, whose goal is to provide first a 2 ((0, ); (Ω)) estimate on the approximate mean chemical potential ,Δ , and then a non-weighted 2 ( ) estimates on the chemical potentials. Proof. Let ≥ 1 and = | ∈ ℰ int , then thanks to (3.20), either 1, ≥ ⋆ 2 or 2, ≥ ⋆ 2 . Let us assume that 1, ≥ ⋆ 2 , the other case being similar. We can also assume without loss of generality that 2, ≥ 2, ≥ 2, . Then the triangle inequality ensures that Using relation (2.6), the second term of the right-hand side rewrites while since 2, ≥ 2, and 1, ≤ 1 ≤ Therefore, using ( + + ) 2 ≤ 3( 2 + 2 + 2 ), we get that where we have set Using Cauchy-Schwarz inequality, we get that ≤ 12 We deduce from 0 ≤ , ≤ 1, from = 2 and from Lemma 3.6 that Besides, Cauchy-Schwarz inequality yields The first term in the right hand side is bounded uniformly w.r.t. thanks to Corollary 3.4. Reorganizing the second term, one gets that Thanks to assumption (2.11) on the regularity of the mesh, one has ∑︀ ∈ℰ ,int ≤ . Therefore, it follows from Lemma 3.5 that (Ω)) estimate (3.23) with the zero mean condition (2.7) allows to make use of a Poincaré-Sobolev inequality (see for instance [21], [17], Lem. 5.1 or [5]). This provides the following uniform 2 ( ) estimate on the discrete global chemical potential (recall here that Ω ⊂ R 2 ): The definition (2.7) of and the equation (2.6) provide the following relations: As a result of Lemmas 3.2, 3.5 and Estimate (3.28), we recover (3.24).

Existence of a discrete solution
We are now in position to finish the proof of Theorem 2.3 by showing the existence of (at least) one discrete solution to the scheme (2.3)-(2.7). Proof. The proof relies on a topological degree argument [12,24]. Our goal is to pass continuously from a linear problem for which the existence and uniqueness of the solution is known to the nonlinear system given by our scheme. Since the construction of such an homotopy (which is parametrized by ∈ [0, 1]) is non-trivial, we give here a description of it, as well as of the key estimates that allow us to use this machinery.
from which we infer that 0 1, = 0 -and thus 0 2, = 0 because of (3.35b), the right-hand side having been set to 0 -for all ∈ and that 0 , does not depend on . Moreover, we deduce from (3.35c) with zero right-hand side that 0 1, = 0 2, for all ∈ . Finally, (3.35d) shows that the discrete chemical potential 0 , are all equal to 0, hence Ker(L) is trivial. Therefore, L has maximal rank and its span can be characterized as the orthogonal of one non-zero vector generating Ker(L ). Denote by 1 = (1, . . . , 1) ∈ R # , 0 = (0, . . . , 0) ∈ R # , and by = ( ) ∈ ∈ R # , then one readily checks that where we have set Elementary convexity inequalities yield and, setting On the other hand, using the boundedness of −1 , between 0 and 1, one gets that while the boundedness of Ψ , yields where depends only on , Ω, ‖Ψ ‖ ∞ , , , and ℎ (but not on ), and where we have set with the convention that ( ) = +∞ if / ∈ [0, 1] and = 1. As a consequence of the technical Lemma A.1 stated in appendix, there exists depending only , , , ℎ , ‖Ψ ‖ ∞ and such that ( Besides, performing a discrete integration by parts on the term 1 yields Since is 1-Lipschitz continuous, one has This implies in particular that is bounded independently uniformly w.r.t. , hence for some not depending on . We can derive a control on for < 1 from the control of the energy dissipation D in (3.39), but this control degenerates as tends to 1. To bypass this difficulty, we multiply (3.30) by (︀ , )︀ . Since , has been designed so that we can mimic the proof of Lemma 3.5 in order to get the existence of not depending on such that Thanks to (3.40), we can reproduce the proof of Lemma 3.7 to claim that for some ⋆ > 0 not depending on . This provides a uniform in discrete estimate on (︀ )︀ ∈ and finally the existence of some 9 not depending on following the path of Lemma 3.8 such that (3.42) Then the topological degree corresponding to system (3.30)-(3.34) on the compact set

Convergence of the scheme
The goal of this section is to prove Theorem 2.4, i.e., that ( ,Δ , ,Δ ) tends to a weak solution ( , ) of (1.1)-(1.6) in a suitable topology as ℎ and ∆ tend to 0 provided the mesh remains sufficiently regular. Consider a sequence of regular meshes (︀ , ℰ , ( ) ∈ )︀ ≥1 such that (2.11)-(2.14) hold for some uniform , ℓ ⋆ , ⋆ and ⋆ w.r.t. , and such that ℎ tends to 0 as tends to +∞, and a sequence of times steps (∆ ) ≥1 with ∆ = / with tending to +∞ with . Then the a priori estimates derived in Section 3.1 are satisfied uniformly provided (3.20) holds, as it is the case if the inverse CFL condition (3.18) is fulfilled.
The first lemma gathers some first consequences of the a priori estimates stated in Section 3.1.
Lemma 4.1. There exist ∈ ∞ ( ; [0, 1]) with 1 + 2 = 1 and ∈ 2 ( ), ∈ {1, 2} such that, up to a subsequence, a.e. in and in the ∞ ( ) weak-⋆ sense, Moreover, ∫︀ Ω ( , ) d = 0 for a.e. ≥ 0, where = 1 1 + 2 2 . Proof. Because of Lemma 3.2, the approximate solutions , ,Δ remain bounded a.e. in between 0 and 1. Therefore, there exists ∈ ∞ ( ; [0, 1]) such that, up to a subsequence, , ,Δ tends to in the ∞ ( )weak-⋆ sense. This is enough to pass in the limit in the relation 1, ,Δ + 2, ,Δ = 1 which directly follows from (2.6). On the other hand, it follows from Lemma 3.8 that the sequences ( , ,Δ ) ≥1 are uniformly bounded in 2 ( ), hence the weak convergence in 2 ( ) towards some . The almost everywhere convergence of , ,Δ towards follows from the discrete Aubin-Lions lemma stated in Appendix B. Finally, given an arbitrary ∈ ∞ (0, ), then multiplying ( ,Δ ] d d = 0. Since , ,Δ converges a.e. towards while remaining uniformly bounded between 0 and 1, it converges also in the strong 2 ( ) sense thanks to Lebesgue dominated convergence theorem. Together with the weak convergence in 2 ( ) of , ,Δ towards , we have enough compactness to pass to the limit → +∞ in the above expression, which gives that This of course implies that ∫︀ Ω ( , ) d = 0 for a.e. ≥ 0. Before going further, we need to introduce some additional material concerning the construction of a strongly consistent approximate gradient based on the SUSHI finite volume scheme [17]. We gather in the following proposition the properties of this approximate gradient to be used in what follows. The super-admissibility of the mesh is crucial at this point. We refer to [17] or to Chapter 13 of [13] for the proofs corresponding to Proposition 4.2.
Then there exists a linear operator ∇ : → ∞ ( ) 2 such that: and all ∈ {1, . . . , }, one has (ii) if the sequence ( ,Δ ) ≥1 is such that ‖ ,Δ ‖ 2 ( ) and ‖∇ ,Δ ‖ 2 ( ) 2 are bounded w.r.t. , then there exists ∈ 2 ((0, ); 1 (Ω)) such that Let us point out that we could have improved the convergence property in (iii) until obtaining the uniform convergence at the price of adding some additional degrees of freedom on the boundary edges. However, the convergence properties stated in Proposition 4.2 are sufficient to prove the convergence of our scheme. Therefore, we avoid the introduction of additional material.
Next lemma focuses on the term corresponding to ∇ . For ≥ 1, we define and the corresponding piecewise constant vector field Proof. Since = 2 and since 0 ≤ , ≤ 1, it results from Lemma 3.6 that Therefore, up to a subsequence, , ,Δ converges weakly in 2 ( ) 2 towards some . Let us identify as − ∇ . To this end, we introduce an arbitrary smooth vector field Φ ∈ ∞ ( ) 2 , and, for all ≥ 1, we denote by and by for almost all ( , ) ∈ . Thanks to the regularity of Φ, it is easy to check that both Φ ,Δ and Φ ,Δ converge uniformly towards Φ as tends to +∞. This implies in particular that On the other hand, , (Φ) can be decomposed into , (Φ) + Let us first focus on (1) , (Φ), which can be controled as follows thanks to Cauchy-Schwarz inequality and the fact that ≤ 2ℎ : The second sum in the right-hand side is uniformly bounded thanks to Lemma 3.6, whereas since | , − , | ≤ | , − , |, Lemma 3.5 and the particular definition (2.4) of , ensure that Since , ≥ ℎ , we finally obtain that Let us now consider (2) , (Φ). To this end, remark first that the definition of Φ implies that As a consequence, since , ,Δ converges weakly towards and , ,Δ converges strongly towards in 2 ( ), we conclude that Thanks to (4.1), the term (3) , (Φ) can be rewritten as The strong convergence of ∇ , ,Δ towards ∇ in 2 ( ) 2 , the weak convergence of , ,Δ towards and the uniform convergence of Φ ,Δ towards Φ yield Introducing the quantities and the corresponding functions , ,Δ in ,Δ , the term (4) , (Φ) rewrites Since , ,Δ is uniformly bounded in 2 ( ), proving that , ,Δ tends to 0 in 2 ( ) is enough to show that Thanks to the regularity of the mesh , and more precisely to (2.13), there holds Using estimate one infers from (2.12) that so that , ,Δ tends to 0 in 2 ( ) and (4.9) holds. Finally, we deduce from (4.5) to (4.9) that = − ∇ in the distributional sense, hence also in 2 ( ) 2 .
We have now all the necessary material to conclude the proof of Theorem 2.4. This is the purpose of our last lemma.
Proof. As a preliminary, let us first show that the functions , ,Δ defined by converges strongly in 2 ( ) towards . Indeed, one has Since , ,Δ converges in 2 ( ) towards as tends to ∞, then so does , ,Δ . Let ∈ ∞ ([0, )×Ω), then denote by = ( , ) for all ∈ and all ∈ {0, . . . , }, ≥ 1. Note that = 0 for all ∈ . Multiplying (2.3) by ∆ −1 and summing over ∈ and ∈ {1, . . . , } leads to , + , + , + , = 0, (4.10) where we have set Classical arguments (see for instance [16]) allow to show that and, since , tends to , that Using Taylor expansions, one shows that (4.13) Therefore, Cauchy-Schwarz inequality together with Lemma 3.6, the regularity of and the 2 ( ) bound of    [30] in (2.3). This scheme degenerates into the upstream mobility scheme proposed in [9] in the deep quench limit , = 0. Almost all our analysis can be adapted to this scheme excepted Lemma 4.5. More precisely, we are not able to prove that the term (1) , (Φ) appearing in the proof of Lemma 4.5 tends to 0 as tends to +∞, which possibly breaks the consistency of the scheme.

Numerical results
In this section, we present different simulations to illustrate the behavior of the finite-volume scheme presented in Section 2.2. To solve this nonlinear system we use a Newton-Raphson based iterative method. More precisely, the unknowns (︀

2,
)︀ ∈ are eliminated thanks to the relation (2.6), so that the nonlinear system to be solved  at each time step involves 3 unknowns 1, , 1, and 2, per cell ∈ . The iterative method stops as soon as the ℓ 2 norm of the Newton increment is smaller than 10 −6 . The updated concentration variables are projected on the set [ , 1 − ] , with = 10 −10 , which is reasonable in view of Lemma 3.2.
In each case the domain Ω is the square (0, 1) 2 . The mesh is made of 23330 conforming triangles. The mesh size is approximately equal to 0.017 and the time step is fixed to ∆ = 10 −4 . We choose as parameters = 0.0002, = 1.45, 1 = 2 = 0.35, = 1 and 1 = 2 = 1. We plot the concentration 1 and we can observe in blue the concentration 1 = 0, in red 1 = 1 and in white 1 = 0.5.
First we consider the spinodal decomposition test case. The initial saturation 0 1 is a random initial concentration with a fluctuation, that is 0 1 ( ) = 0.5 + ( ) where ≪ 1 is a small random perturbation. We compare the case without any external potential, that is Ψ 1 = Ψ 2 = 0, in Figure 2 with the case where the external potential are given by Ψ ( ) = − · where the gravity is = −0.98 and the densities 1 = 5 and 2 = 1 in Figure 3. Note that in both cases we have exactly the same initial data. We want to observe the influence of the gravity on the phase separation dynamics.
At the very beginning (see Figs. 2a and 3a), as the state 1 = 0.5 is slightly disturbed, the two pure phases   Now we keep exactly the same data but we change the initial condition by favouring the pure phase 1 = 0 and choosing 0 1 ( ) = 0.3 + ( ). First of all we can see in Figures 4a and 5a that, as expected, the pure phase 1 = 1, and thus, the phase separation dynamics appears later than in Figures 2a and 3a. Moreover, since the phase 1 = 0 is preponderant, a collection of circular droplets of the pure phase 1 = 1 can be observed over a long period of time (see Figs. 4c, 4d and 5c, 5d). Our observations on Figures 2 and 3 concerning the influence of the external potentials are still valid with this choice of initial profiles.
We consider now a second test case. The initial concentration is a cross in the middle of the domain presented in Figure 6a. Here again we start with the case without external potentials. We know that in this case the Cahn-Hilliard model preserves the volume while minimizing the perimeter and thus as seen in Figure 6 the cross evolves into a circle. Now, we want to observe the influence of the gravity when we add the external potentials Ψ ( ) = − · . Since 1 = 5 > 1 = 2 , as one might expect, we observe in Figure 7 that the cross, that is the pure phase 1 = 1, is drawn down. Thus, although the volume is still preserved, the final state is no longer a circle but a strip at the bottom of the domain. Let us establish the following lemma, which is used in the proof of Proposition 3.9.