On the viscous Cahn-Hilliard-Navier-Stokes equations with dynamic boundary conditions

In the present article we study the viscous Cahn-Hilliard-Navier-Stokes model, endowed with dynamic boundary conditions, from the theoretical and numerical point of view. We start by deducing results on the existence, uniqueness and regularity of the solutions for the continuous problem. Then we propose a space semi-discrete finite element approximation of the model and we study the convergence of the approximate scheme. We also prove the stability and convergence of a fully discretized scheme, obtained using the semi-implicit Euler scheme applied to the space semi-discretization proposed previously. Numerical simulations are also presented to illustrate the theoretical results.

1. Setting of the problem. In this paper we consider the coupling between the incompressible Navier-Stokes equations with the viscous Cahn-Hilliard equations and we endow the Cahn-Hilliard equations with dynamic boundary conditions. This coupling of equations models the motion of isothermal mixture of two confined immiscible and incompressible with comparable (equal) densities and viscosities; when the two fluids have comparable constant densities the temperature differences are neglectable and the diffusive interface between the two phases has a small, non-zero thickness. The model reads as follows: in Ω T , ϕ t + u · ∇ϕ − γ∆w = 0 in Ω T , w = −∆ϕ + f (ϕ) + εϕ t , in Ω T , u = 0, ∂ n w = 0 on Γ × (0, T ), with ε > 0 a viscosity parameter, δ > 0, λ s > 0, γ > 0 the elastic relaxation time, λ > 0 the surface tension and Ω a bi-dimensionnal domain with smooth boundary Γ. Here u and p respectively denote the average velocity and the pressure of the fluid mixture and the scalar ϕ is an order parameter indicating the phases of the fluid. The two fluids are supposed to have the same density, taken here equal to 1, and the same viscosity, denoted by ν. In (1), ∆ Γ is the Laplace-Beltrami operator on Γ, n and ∂ n are respectively the outer unit normal and the associated normal derivative on Γ. The evolution boundary value problem (1) is completed by the initial condition u(0) = u 0 and ϕ(0) = ϕ 0 .
The functions f and g s are polynomials of degree p and q respectively, with p ≥ 3 odd, q ≥ 1 odd and the dominant coefficients positive. Typical choices are We will also consider F and G s to be two antiderivatives of the functions f and g s . The Cahn-Hilliard-Navier-Stokes (CHNS) model with classical (non-dynamical) boundary conditions was highly studied from the theoretical and mathematical point of view. The model was first introduced in physical papers as [6], [8] (see also [20,24] where the model was derived using mathematical arguments). The mathematical study of the problem, meaning the existence and uniqueness of solutions for the two and three dimensional cases, as well as the asymptotic stability of the model was established in [3]. In [15] the authors studied the asymptotics of the CHNS model proving the existence of a global finite dimensional attractor those fractal dimension is estimated; the convergence of each trajectory to an equilibrium is established. The case of a CHNS model with degenerate mobility for the Cahn-Hilliard equation was also studied in [1]. The study of pullback attractors for Cahn-Hilliard-Navier-Stokes models can be found in [2,26].
For the numerical studies, there is also a lot of litterature where different kind of approximate models are considered. We refer the interested reader to [11,7,12,13,21,22], where different finite element type of approximations for interface models are considered.
The novelty in this paper is that we consider the case of a coupling between the viscous Cahn-Hilliard equations endowed with dynamic boundary conditions and the Navier-Stokes equations. To the best of our knowledge, only one work concerns the same subject, see [16], where the authors consider the theoretical study of the problem, meaning the well posedness of the model. A study of bounded liquid-gas flows with different dynamic boundary conditions can be found in [9].
The paper is structured as follows: the first two sections address the theoretical study of the model while the rest of the paper proposes discrete approximations, first obtained by the discretisation of the space variable using a finite element approximation and then by introducing a fully discrete model. Thus, in Section 2 we derive a priori estimates on which the existence, uniqueness and regularity of the solutions are proved. In Section 3 the existence of a weak solution is rigorously derived. In Section 4 we propose a space semi-discrete finite element approximation for the continuous problem and in Section 5 we derive error estimates and we prove the convergence of the sequence of approximate solutions to the exact solution. Section 6 is devoted to the study of a fully discretized scheme, obtained using the backward Euler method for the discretization of the time derivative in the space semi-discrete scheme studied previously. The existence of a unique solution for the fully discretized model, the conditional stability and the convergence of the method are established. In the end we also present some numerical simulations that illustrate the results presented above.
Using the fact that F (ϕ) ≥ −c 1 and G s (ϕ) ≥ −c 2 , where c 1 , c 2 are positive, absolute constants, the first three lines of Lemma 2.1 are direct consequences of (8) and (9). We also infer that ∇w ∈ L 2 (0, T ; L 2 (Ω)). Thus, we only need an estimate on w in order to be able to conclude that the H 1 − norm of w remains bounded. From (3) 3 , taking χ = 1 and recalling that ϕ t = 0, we obtain: Using (10), we conclude that w belongs to L 2 (0, T ) and thus w ∈ L 2 (0, T, H 1 (Ω)). It remains to detail the estimate for p. To this aim, we integrate (4) with respect to t ∈ [0, τ ], τ ≤ T . This yields thanks to (8) and (9). Recalling the Inf-Sup condition (see [18]): we get This completes the proof of the lemma.
We take v = u in (11) 1 . Using the fact that Now we take ψ = w in (11) 2 and χ = ϕ t in (11) 3 and sum the resulting equations. This yields to (13) Using the fact that f and g s are polynomials, we can write: with the constant c depending on 1/ε. We also have: Moreover, since ∇ · u i = 0 for i = 1, 2, we have: and we bound the terms from the right hand side of (14) as follows: Similarly, we compute: Hence, setting we infer from (12), (13) and the estimates above that the following estimate holds: where the function h belongs to L 1 (0, T ). Thus, we deduce from the Gronwall Lemma that: H(t) ≤ cH(0) = 0 ∀t ≥ 0 (since ϕ(0) = 0 and u(0) = 0).
In particular, we infer that, for every t ≥ 0, which implies the uniqueness of a weak solution.

More regularity results.
In what follows we prove the following higher order regularity results for the solution of (1): Lemma 2.2. We assume that u 0 ∈ (H 1 0 (Ω)) 2 , ϕ 0 ∈ H 1 (Ω, Γ). Then the solution to (1) satisfies, for all T > 0: We first prove the H 2 -regularity for ϕ. To this aim, we rewrite the fourth and sixth equations of (1) as with the right hand side of these equations belonging respectively to L 2 (0, T ; L 2 (Ω)) and L 2 (0, T ; L 2 (Γ)). Then, we apply the Maximum Principle introduced in [25] and obtain: . Now, having the L ∞ -estimates for ϕ and ϕ |Γ , we can interpret the nonlinearities f (ϕ) and g s (ϕ) as external forces as well and apply the H 2 -regularity to the linear elliptic system and obtain that ϕ is bounded in L 2 (0, T ; H 2 (Ω, Γ)) (for more details, see [25]). Now we focus on the regularity for w. We write (1) 3 as: Noticing that u · ∇ϕ L 3/2 (Ω) ≤ |u| Ω |∇ϕ| L 6 (Ω) , we obtain in view of Lemma 2.1 and of the above proved regularity for ϕ. Using the fact that ϕ t + u · ∇ϕ and w are in L 2 (0, T, L 3/2 (Ω)), we apply the H 2 -regularity for an elliptic linear system and we obtain w ∈ L 2 (0, T ; W 2, 3 2 (Ω)). Now we focus on the regularity for u. We take v = Au in (3) 1 , where A is the Stokes operator (see [28]). We obtain: We can write: and, using the fact that W 1, 3 2 (Ω) → L 6 (Ω), Gathering all the estimates, we obtain the following energy estimate: , we can apply the Gronwall lemma to (18) and we obtain that u ∈ . Considering again problem (16), we now have ϕ t + u · ∇ϕ in L 2 (0, T ; L 2 (Ω)). Then, applying one more time the H 2 -regularity, this implies w ∈ L 2 (0, T ; H 2 (Ω)).
Proof. We differentiate (1) with respect to t and find that (u t , p t , ϕ t , w t ) satisfies the following system: We multiply the first equation of (19) by u t and integrate over Ω. It yields: We first notice from (5) that ϕ t = 0. Multiplying (19) 3 by w t , (19) 4 by ϕ tt , integrating over Ω and summing the resulting equations, we obtain: Now we notice that (since ∇ · u t = 0 and u t = 0 on Γ) Therefore, we sum (20) to equation (21) multiplied by λ, and use the following bounds: with and Applying Lemma 2.1, we infer that Λ defined in (24) belongs to L 1 (0, T ). Hence, we can apply Gronwall Lemma to get the announced regularity for u t , ϕ t and ϕ tt . Since u t , u ∈ L 2 (0, T ; (H 1 0 (Ω)) 2 ), then using Lemma 1.2 in [23] and eventually modifying u on a set of zero measure in (0, T ), we can conclude that u ∈ C 0 ([0, T ]; (H 1 0 (Ω)) 2 ).

Remark 2.
We emphasize that the constant C appearing in (24) depends on 1/ε. Hence we will not be able to obtain results for the zero viscosity Navier-Stokes-Cahn-Hilliard problem (i.e. ε = 0).
3. Existence of solutions. We consider the following problem, in order to obtain the Galerkin approximation for u: This problem defines the operator A 1 such that is a compact, self-adjoint operator. Hence there exists (f j ) j ⊂ W 0 a sequence of orthonormal eigenfunctions such that A 1 f j = λ j f j , λ j > 0; the family (f j ) j spans (L 2 (Ω)) 2 and satisfies, for all i, k ∈ N: We also introduce the operator N which is the inverse of the Laplacian operator (−∆), where −∆ is endowed with Neumann boundary conditions, imposing zero average over the domain Ω. We take {e n } n orthonormal in L 2 (Ω), e 1 = const (β 1 = 0) and N e i = 1 βi e i , β i > 0, for all i ≥ 2. Note that e i = 0 ∀i ≥ 2. We introduce the following finite dimensional approximate problem: . . , f n }, ψ, χ ∈ V n = Span{e 1 , e 2 , . . . , e n }, where P n is the projection operator from (L 2 (Ω)) 2 into W n . Now we take v = f j , j = 1 . . . n and χ = ψ = e j , j = 1 . . . n. We obtain: We define the matrices: and the functions We also introduce the vector functions defined by: (H) j = (P n h, f j ) Ω , ∀ j = 1, . . . , n.' Then, problem (26) becomes the following system of nonlinear ordinary differential equations: Combining the last two equations in (27), we obtain: In order to be able to apply the Cauchy-Lipschitz theorem to (27) 1 and (28), we need to prove that the matrix N = I M + εγM Ω + γM Ω A Γ is invertible. Since e 1 = const, we emphasize that M Ω has all elements from the first line and the first column equal to zero. As a consequence, M Ω A Γ has its first line equal to zero, and we are lead to prove thatṄ is invertible, whereṄ is the (M − 1) × (M − 1) matrix obtained, deleting the first line and the first column of N . The matrixṀ Ω is diagonal, with all terms on the diagonal strictly positive. In particular,Ṁ Ω is positive definite, and so is the matrixṀ −1 Ω + εγI M −1 . We also notice that: N = I M −1 + εγṀ Ω + γṀ ΩȦΓ =Ṁ Ω (Ṁ −1 Ω + εγI M −1 + γȦ Γ ). SinceȦ Γ is positive semi definite, we conclude thatṄ is invertible. Thus we can solve the system in (Φ, U ). By the Cauchy-Lipschitz theorem, problem (25) has a unique maximal solution (u n , ϕ n , w n ) in C 1 ([0, T + ); W n × V n × V n ).

Finite element approximation.
In this section and all the following sections, the domain Ω will be considered to be a 2d slab, meaning Ω = (0, L 1 )×(R/L 2 Z) with smooth boundary Γ = {0, L 1 } × (R/L 2 Z). The purpose of this section is to propose a finite element approximation for the exact solution of (1), where all the unknowns in (1) are searched to be periodic in the x 2 direction, while the Dirichlet boundary condition for u and the dynamic boundary condition for (ϕ, w) are imposed on x 1 = 0 and x 1 = L 1 .
Starting with this section, the typical function spaces are: where by H m p we understand the functions that belong to H m (Ω) and which are periodic in the x 2 -direction, while H m 0,p (Ω) = {v ∈ H m p (Ω); v = 0 at x 1 = 0 and x 1 = L 1 }. We remark here that all the results proved previously still remain valid when we consider the domain with periodic boundary conditions in one direction.
Let J h be a quasi-uniform partitioning of Ω into disjoint open simplices K. We set h K = diam (K) and h = max K∈J h h K . The mesh J h is assumed to be weakly acute, that is for any pair of adjacent triangles the sum of the opposite angles relative to the common side does not exceed π. We also consider J h/2 to be the mesh obtained from J h by refining each simplex K into four similar triangles obtained joining the midpoints of each edge of K. We The triangulation J h ofΩ induces in a natural way a triangulation Γ h of Γ. Note that for every v h ∈ V h , the restriction v h | Γ on the boundary is a P 1 -finite element on the one dimensional domain Γ. Thus, V h is also used for the discretisation of H 1 p (Ω, Γ) (see [4] for more details). We now define the projections Q 1 h : (L 2 (Ω)) 2 → (V h/2 ) 2 , such that: and Q 2 h : L 2 (Ω) → V h , such that: We introduce the elliptic projections of ϕ and w on V h : for all φ ∈ V h . Then, the following result holds (see [5] for more details): Lemma 4.1. For all ϕ ∈ H 2 p (Ω, Γ) and for all w ∈ H 2 p (Ω), the functions (φ h ,w h ) ∈ V h × V h defined by (34), satisfy: and where | · | m,Ω and | · | m,Γ are the seminorms associated with respectively the H m p (Ω) and H m p (Γ) norms.
We now introduce the following finite element spaces: and the elliptic projection of u on W h 0 which is defined by: We emphasize that the space W h 0 is continuously embedded in (H 1 0,p (Ω)) 2 , but not in W 0,p . The following result is classical: For all v ∈ W 0,p , the projectionṽ h defined by (37) satisfies: The semi-discrete scheme reads (we take λ = 1): for all v ∈ W h 0 and ψ, χ ∈ V h . Taking ψ = 1 in (38) 2 , we obtain: We thus obtain the mass conservation for the approximation ϕ h : (ϕ h , 1) Ω = const ∀t ≥ 0.
We take v = u h in (38) 1 , ψ = w h in (38) 2 and χ = ϕ h t in (38) 3 and we get: with J as in (7). This gives the following lemma: The semi-discrete scheme satisfies the following stability estimate: The existence of the solution for the discrete problem (38) follows from Lemma 4.3, using the theory of ordinary differential equations.

5.
Error estimates for the semi-discrete scheme. In order to estimate the error between the exact solution of the variational problem (3) and the approximate solution given by (38), we use a splitting method introduced by [10] (see also [17], [5] for the application of the splitting method in different contexts). We write: where (ũ h ,w h ,φ h ) are the elliptic projections of (u, w, ϕ) defined by (37), (34).
We also define the discrete inverse Laplacian T h : (39) Note that (39) still holds for χ h ∈ V h since f is a function of zero average.
We also define the discrete negative seminorm: . Taking into account Lemma 4.1 and 4.2, we already have estimates on ρ u , ρ ϕ , ρ w .
It remains to estimate θ u , θ ϕ , θ w ; in order to do so we write the equations that they satisfy: (4) and (38), we obtain the following equation for θ u : In order to deduce a priori estimates on θ u , we take v = θ u in (40) and obtain: (41) We need to estimate the terms from the right hand side of (41). We have: where we used · −1,Ω to denote the norm on H −1 (Ω) = (H 1 0 (Ω)) . We also used that ρ u t −1,Ω ≤ ch |u t | Ω . We also know that: h p ∈ V h and θ u ∈ W h 0 , and we can write: For the nonlinear terms in b, we proceed as follows: Thanks to the orthogonality property b(u, v, v) = 0, we obtain: The terms from the right hand side of the latter inequality are estimated as follows: since, using the fact that u ∈ C 0 ([0, T ], (H 2 p (Ω)) 2 ), the term |∇ũ h | Ω can be estimated as follows: We also estimate: and for each term we proceed as follows: |∇θ w | 2 Ω + c|θ u | 2 Ω , and similarly, |(ϕ∇ρ w , θ u ) Ω | ≤ |∇ρ w | 2 Ω + c|θ u | 2 Ω , where we used that ϕ ∈ L ∞ (0, T ; H 2 (Ω)).
Proof. The proof follows exactly the same lines as Lemma 4.4 in [21] and will be omitted.
For simplicity, the problem is solved with a slightly different discretization (semiimplicit scheme for time discretization, combined with a finite element for the space discretization: P 1 -finite element for ϕ, w and p, and P 2 for u).
In the figures below, the negative and the positive values of the solution ϕ are respectively displayed in yellow and red.
In our first test, we take u 0 and ϕ 0 as in Figure 1. Then, Figures 2 and 3 illustrate different values of ϕ, corresponding to the evolution of the order parameter with respect to time: t = 0.5, t = 1, t = 2.5 (in Figure 2), t = 3, t = 4 and t = 7.5 (in Figure 3). In Figure 4, we show the solution u and p corresponding to t = 3.
For our second test, we take ϕ 0 randomly distributed between −1 and 1, and u 0 and h are again defined as in (73). We compare the solution ϕ at t = 1 ( Figure  5) and t = 3 ( Figure 6), but with different kind of boundary conditions on the horizontal sides of Ω, namely: dynamic boundary conditions withg s = 2v + 1 (left column), dynamic boundary conditions withg s = 2v − 1 (middle column) and Neumann boundary conditions (right column). Figure 7 is obtained taking again ϕ 0 randomly distributed between −1 and 1, and u 0 and h defined as in (73), and we choose the surface potentialg s = 2v + 1. Then we compare the solutions ϕ at the same time t = 3, but with three different values for the viscosity coefficient: ε = 0.5, ε = 0.1 and ε = 0.01.
We also remark that checking the first two columns of Figure 5 and 6, corresponding to the case of dynamic boundary conditions, the contact angle to the boundary is not π/2, contrary to the Neumann case. Nevertheless, in the present work we did not control the contact angle; this question will be addressed elsewhere, eventually with different dynamic boundary conditions.