A linear, second-order, energy stable, fully adaptive finite-element method for phase-field modeling of wetting phenomena

TL;DR

This paper introduces a new finite-element numerical method for phase-field modeling of wetting phenomena that is energy stable, adaptive, and preserves mass, demonstrated through various realistic tests.

Contribution

It presents a novel, fully adaptive, second-order finite-element method for the Cahn-Hilliard equation with wetting boundary conditions, ensuring energy stability and mass conservation.

Findings

Method is mass-conservative and energy-stable

Accurately models wetting on heterogeneous substrates

Effective in complex 2D and 3D scenarios

Abstract

We propose a new numerical method to solve the Cahn-Hilliard equation coupled with non-linear wetting boundary conditions. We show that the method is mass-conservative and that the discrete solution satisfies a discrete energy law similar to the one satisfied by the exact solution. We perform several tests inspired by realistic situations to verify the accuracy and performance of the method: wetting of a chemically heterogeneous substrate in three dimensions, wetting-driven nucleation in a complex two-dimensional domain and three-dimensional diffusion through a porous medium.

Tables

Table 1. Table 1 : CPU times (hh:mm:ss) using an Intel i7-3770 processor for the simulations presented in Section   4.3 (two droplets on a substrate), with or without time-step adaptation. The method OD2-W was used for the simulations with an adaptive time step, and the method OD1-W was used for the simulations with a fixed time step. In both cases, an adaptive mesh was used, with the parameter h min subscript ℎ h_{\min} equal to ε / 10 = 0.001 𝜀 10 0.001 \varepsilon/10=0.001 .

Contact angle	Adaptive time step	Fixed time step
$\pi /4$	44:15:17	130:16:16
$\pi /2$	21:38:38	128:35:20
$3\pi /4$	31:12:18	122:33:28

Table 2. Table 2 : CPU times (hh:mm:ss) using an Intel i7-3770 processor for the simulations presented in Section   4.4.3 (nucleation in a domain with complex boundaries), with or without time-step adaptation. The method OD2-W was used for the simulations with an adaptive time step, and the method OD1-W was used for the simulations with a fixed time step. In all cases, we used a fixed mesh with mesh size h = 0.01 ℎ 0.01 h=0.01 (the size of the domain is roughly 5 by 5) and P ​ 1 𝑃 1 P1 elements, leading to 181587 unknowns. The parameter ε 𝜀 \varepsilon was set to 0.02 0.02 0.02 .

Contact angle	Adaptive time step	Fixed time step
$\pi /4$	36:13:15	65:48:44
$\pi /2$	32:58:01	65:24:10
$3\pi /4$	36:43:13	67:00:50

Equations

$$E(\varphi )$$

$$:={\int }_{\Omega }\left(\frac{1}{\epsilon }{F}_m(\varphi )+\epsilon \frac{|\nabla \varphi {|}^2}{2}\right)\mathrm{d}\Omega +{\int }_{\partial \Omega }{F}_w(\varphi )\mathrm{d}\sigma ,$$

$$F_w(\varphi )=\frac{\sqrt{2}}{2}\cos \theta (\mathbf{x})\left(\frac{{\varphi }^3}{3}-\varphi \right),$$

$$\displaystyle\mathchoice{\frac{\mathrm{d}{}}{\mathrm{d}\alpha}}{\partial_{\alpha}}{}{}E(\phi+\alpha\,\psi)\big{|}_{\alpha=0}$$

$$\mu :=\frac{\delta E}{\delta \varphi }=\frac{1}{\epsilon }{f}_m(\varphi )-\epsilon \Delta \varphi ,$$

$$\epsilon \nabla \varphi \cdot \mathbf{n}=-f_w(\varphi )=\frac{\sqrt{2}}{2}\cos \theta (\mathbf{x})(1-{\varphi }^2).$$

$$\frac{\partial \varphi }{\partial t}=\nabla \cdot (b(\mathbf{x})\nabla \mu ),$$

$$\frac{\mathrm{d}}{\mathrm{d}t}M(\varphi ):=\frac{\mathrm{d}}{\mathrm{d}t}{\int }_{\Omega }\varphi \mathrm{d}\Omega ={\int }_{\partial \Omega }b\nabla \mu \cdot \mathbf{n}\mathrm{d}\sigma ,$$

$$\frac{\partial \varphi }{\partial t}$$

$$\mu$$

$$\epsilon \nabla \varphi \cdot \mathbf{n}$$

$$b\nabla \mu \cdot \mathbf{n}$$

$$\frac{\mathrm{d}}{\mathrm{d}t}E(\varphi (t))=-\Vert \sqrt{b}\nabla \mu {\Vert }_{L^2(\Omega )}^2+{\int }_{\partial \Omega }\dot{m}\mu \mathrm{d}\sigma .$$

$$\varphi (\mathbf{x})=\tanh \left(\frac{\mathbf{x}\cdot \mathbf{u}}{\sqrt{2}\epsilon }\right),\text{where }\mathbf{u}=(\pm \sin \theta ,\cos \theta {)}^T$$

$${\Vert {\varphi }^3\Vert }_{L^1(\partial \Omega )}$$

$$=C({\Vert {\varphi }^3\Vert }_{L^1(\Omega )}+3{\Vert {\varphi }^2\nabla \varphi \Vert }_{L^1(\Omega )})$$

$$\le C({\Vert {\varphi }^3\Vert }_{L^1(\Omega )}+\frac{3}{2\alpha }{\Vert {\varphi }^4\Vert }_{L^1(\Omega )}+\frac{3\alpha }{2}{\Vert \nabla \varphi \Vert }_{L^2(\Omega )}^2)\forall \alpha >0,$$

$$F_w^{\ast }(\varphi )=\frac{\sqrt{2}}{2}\cos \theta (\mathbf{x})\times \{\begin{array}{ll}\left(\frac{2}{3}-(\varphi +1{)}^2\right)& \text{ if }\varphi <-1;\\ \left(\frac{{\varphi }^3-3\varphi }{3}\right)& \text{ if }\varphi \in [-1,1];\\ \left(-\frac{2}{3}+(\varphi -1{)}^2\right)& \text{ if }\varphi >1,\end{array}$$

$$\varphi \in L^{\infty }(0,T;H^1(\Omega )),\frac{\partial \varphi }{\partial t}\in L^2(0,T;(H^1(\Omega ){)}^{\prime }),\mu \in L^2(0,T;H^1(\Omega )),$$

$$⟨{\partial }_t\varphi ,\psi ⟩+\left(b\nabla \mu ,\nabla \psi \right)={\left(\dot{m},\psi \right)}_{\partial \Omega }$$

$$\left(\mu ,\nu \right)=\epsilon \left(\nabla \varphi ,\nabla \nu \right)+\frac{1}{\epsilon }\left(f_m(\varphi ),\nu \right)+{\left(f_w^{\ast }(\varphi ),\nu \right)}_{\partial \Omega }$$

$$\left({\delta }_t{\varphi }^{n+1},\psi \right)$$

$$\left({\mu }^{n+\frac{1}{\alpha }},\nu \right)$$

$$+\frac{1}{\epsilon }\left({\hat{f}}_m({\varphi }^n,{\varphi }^{n+1}),\nu \right)+{\left({\hat{f}}_w({\varphi }^n,{\varphi }^{n+1}),\nu \right)}_{\partial \Omega }$$

$${\varphi }^{n+\frac{1}{\alpha }+\beta }:=\left(1-\frac{1}{\alpha }-\beta \right){\varphi }^n+\left(\frac{1}{\alpha }+\beta \right){\varphi }^{n+1},$$

$$\frac{{\varphi }^{n+1}-{\varphi }^n}{\Delta t}.$$

$$\frac{E({\varphi }^{n+1})-E({\varphi }^n)}{\Delta t}=-\Vert \sqrt{b}\nabla {\mu }^{n+\frac{1}{\alpha }}{\Vert }_{L^2(\Omega )}^2-ND({\varphi }^n,{\varphi }^{n+1}),n=0,1,\dots ,$$

$$ND({\varphi }^n,{\varphi }^{n+1})=\epsilon N{D}_{\text{philic}}({\varphi }^n,{\varphi }^{n+1})+\frac{1}{\epsilon }N{D}_{\text{phobic}}({\varphi }^n,{\varphi }^{n+1})+N{D}_{\text{wall}}({\varphi }^n,{\varphi }^{n+1}),$$

$$N{D}_{\text{philic}}({\varphi }^n,{\varphi }^{n+1})$$

$$N{D}_{\text{phobic}}({\varphi }^n,{\varphi }^{n+1})$$

$$N{D}_{\text{wall}}({\varphi }^n,{\varphi }^{n+1})$$

$$N{D}_{\text{phobic}}({\varphi }^n,{\varphi }^{n+1})$$

$$\displaystyle-\frac{1}{6}\,\Delta t^{2}\,f_{m}^{\prime\prime}(\phi^{n})\,(\delta_{t}\phi^{n+1})^{2}-\frac{1}{24}\,\Delta t^{3}\,f_{m}^{\prime\prime\prime}(\phi^{n})\,(\delta_{t}\phi^{n+1})^{3}\Big{)}\,\mathrm{d}\Omega,$$

$$N{D}_{\text{wall}}({\varphi }^n,{\varphi }^{n+1})$$

$$-\frac{1}{2\Delta t}{\int }_{{\varphi }^n}^{{\varphi }^{n+1}}{f}_w^{\prime \prime }(\varphi )(\varphi -{\varphi }^n{)}^2\mathrm{d}\varphi \mathrm{d}\sigma .$$

$${\hat{f}}_m({\varphi }^n,{\varphi }^{n+1})$$

$$=\left[\frac{3}{2}({\varphi }^n{)}^2{\varphi }^{n+1}-\frac{1}{2}({\varphi }^n{)}^3-\frac{{\varphi }^n+{\varphi }^{n+1}}{2}\right],$$

$${\hat{f}}_w({\varphi }^n,{\varphi }^{n+1})$$

Full text

Keywords: Wetting, diffuse interface theory, finite element method, Cahn-Hilliard equation, adaptive time step.

1 Department of Chemical Engineering, Imperial College London, London SW7 2AZ, United Kingdom

2 Department of Mathematics, Imperial College London, London SW7 2AZ, United Kingdom

3 School of Mathematics and Statistics, The Open University, Milton Keynes MK7 6AA, United Kingdom

1 Introduction

Capillarity and wetting phenomena, driven primarily by interfacial forces, are ubiquitous in a wide spectrum of natural phenomena and technological applications. Examples range from the wetting of plant leaves by rainwater and insects walking on water to coating processes, inkjet printing, oil recovery and microfluidic devices; for reviews, see e.g. [17, 5]. From an historical point of view, two of the concepts essential to the understanding of capillarity and wetting were introduced and studied already in 1805: these are the Laplace pressure [40] and the Young–Dupré contact angle [67]. Later, following the work of Plateau on soap films [47], Poincaré [48] linked interfacial phenomena with the theory of minimal surfaces.

Wetting phenomena typically involve a fluid-fluid interface advancing or receding on a solid substrate and a contact line formed at the intersection between the interface and the substrate. The wetting properties of the substrate determine to a large extent the behavior of the fluids in the contact-line region, and in particular the contact angle at the three-phase conjunction, defined as the angle between the fluid-fluid interface and the tangent plane at the substrate. At equilibrium, this is precisely the Young–Dupré angle. When one of the two fluids moves against the other, the contact angle becomes a dynamic quantity, and when the problem is formulated in the framework of conventional hydrodynamics, the contact line motion relative to the solid boundary results in the notorious stress singularity there, as first noted in the pioneering studies by Moffat [43] and Huh and Scriven [34]. Since then there have been numerous analyses and discussions of the singularity over the years, see e.g. [61, 32, 15] and also recent studies in [58, 57] (with the latter one revisiting the classical Cox–Hocking matched asymptotic analysis and providing a correction to it). Recently, it was shown [45] that mesoscopic approaches such as dynamic density functional theory (DDFT) can fix this singularity behavior.

A popular model for interface dynamics is the Cahn–Hilliard (CH) equation [10, 12], which belongs to the class of phase-field and diffuse interface models. Originally proposed to model spinodal decomposition, the mechanism by which a binary mixture can separate to form two coexisting phases due to, e.g., a change of temperature [12], it has been used in a wide spectrum of different contexts since, such as solidification phenomena [14] and Saffman–Taylor instabilities in Hele–Shaw flows [31]. To account for wetting phenomena and contact lines on solid boundaries, the CH equation can be coupled to a wall boundary condition [11]. Such CH model has been employed successfully in various situations, including microfluidic devices [18, 19, 49, 65], flow in porous media [4], rheological systems [8], and patterning of thin polymer films [38]. Other potential applications include micro-separators [50], fuel cells [3] and CPU chip cooling based on electro-wetting [44]. Many of these applications are characterized by the presence of chemically heterogeneous substrates and/or complex geometries, which make their numerical simulation challenging.

The form of the wetting boundary condition is dictated by the form of the wall free energy. For liquid-gas problems linear forms have been adopted, e.g. in the pioneering study by Seppecher [51] and in [9, 66]. But a cubic is the lowest-order polynomial required such that the wall free energy can be minimized for the bulk densities, and it prevents the formation of boundary layers on the wall. Cubic forms have been adopted for binary fluid problems, e.g. in [35, 69], but also for liquid-gas ones, see [55, 56]. The latter studies, in particular, showed asymptotically that a CH model can alleviate the contact line discontinuity without any additional physics (and at the same time completing but also correcting Seppecher’s work). The detailed asymptotic analysis of the unification of binary-fluid CH models can be found in [54].

Various approaches have been proposed in the literature for the numerical solution of the CH equation. Because of the high order of the equation and its multiscale features (scale separation between interface size and the characteristic length), most existing time-stepping schemes are implicit or semi-implicit. Several of these schemes aim to satisfy discrete mass and energy laws in agreement with the underlying continuum model. Discretization in space can be achieved using finite-difference methods [25, 39], finite element methods [2, 21, 60], spectral methods [64], or pseudospectral methods [46]. In addition, the computation time can be reduced by applying adaptive mesh refinement [68, 59] and time-step adaptation [28].

Among the several linear schemes for the CH equation with homogeneous Neumann boundary conditions introduced in [27], the authors showed by means of numerical experiments that their second-order optimal dissipation scheme, referred to as OD2, is the most accurate and the one introducing the least numerical dissipation. In this work, we outline a numerical scheme that extends and appropriately generalizes OD2 as follows: (a) it includes a non-linear wetting boundary condition; (b) it adopts an efficient energy-based time-step adaptation strategy. In contrast with the time-adaptation scheme introduced in [28], where the time step is adapted to limit numerical dissipation, we base the time-step adaptation directly on the variation of free energy. With this method we are able to solve the CH system efficiently and systematically to capture wetting phenomena in both two- and three-dimensional (2D and 3D, respectively) settings, and in a wide range of situations, including confinement with complex geometry, chemical and topographical heterogeneities, or both.

Like the OD2 scheme on which it was based, the time-stepping scheme we propose is semi-implicit and linear. We show that it is mass-conservative and satisfies a discrete free-energy law with a numerical dissipation term of order 2 in time. Space discretization is achieved using a finite element method, leading to an unsymmetrical sparse linear system to solve at each iteration. In addition to adapting the time step as mentioned above, with the aim of increasing the resolution in time during fast phenomena, we use a mesh refinement strategy to capture interfaces precisely.

To test the efficiency of the proposed numerical scheme we consider several wetting problems as test cases. We first study relaxation towards equilibrium in two situations: the spreading of a sessile droplet and the coalescence of two sessile droplets on a flat, chemically homogeneous substrate. We then consider two-component systems in complex geometries delimited by chemically heterogeneous substrates in both 2D and 3D.

In Section 2, we introduce the CH system and the non-linear wetting boundary condition. In Section 3, we outline our numerical scheme and prove the associated conservation properties. In Section 4, we present the results of several numerical experiments. Conclusions and perspectives for future work are offered in Section 5.

2 Phase-field model for wetting phenomena

Throughout this study, $\Omega\subset^{d}$ denotes a $d$-dimensional domain, $\partial\Omega$ denotes its boundary with outward unit normal vector n, $\Gamma_{S}$ is the solid substrate and $\Gamma_{G}=\partial\Omega\setminus\Gamma_{S}$. The CH system we use to describe the dynamics of two immiscible fluids in contact with a solid substrate, is a free-energy-based model. The starting point is the introduction of a locally conserved field, denoted by $\phi:\Omega\rightarrow$, that plays the role of an order-parameter: two equilibrium values, say $+1$ and $-1$, represent the pure phases, and the interface is conventionally located at the points where $\phi=0$ [10, 12]. We consider systems with a free energy given by

[TABLE]

where the two terms, $E_{m}$ and $E_{w}$, represent the mixing and wall components of the free energy, respectively. Here $F_{m}(\phi)=\frac{1}{4}(\phi^{2}-1)^{2}$ and $F_{w}$ is taken to be a cubic polynomial, following e.g. [55, 56]:

[TABLE]

where $\theta=\theta(\boldsymbol{\mathbf{x}})$ is the equilibrium contact angle, which can depend on the spatial position $\boldsymbol{\mathbf{x}}$. From the expression of the free energy, we calculate that, for a sufficiently smooth function $\psi:\Omega\rightarrow$:

[TABLE]

with $f_{m}=F_{m}^{\prime}$ and $f_{w}=F_{w}^{\prime}$, so the chemical potential is equal to

[TABLE]

and the natural boundary condition associated with the surface energy is

[TABLE]

We assume that the dynamics of the system is governed by the CH equation,

[TABLE]

where $b(\boldsymbol{\mathbf{x}})$ is a mobility parameter, assumed to be uniform hereafter. This leads to the following mass-conservation property:

[TABLE]

so the mass flux at the boundary can be specified using the condition $b\boldsymbol{\mathbf{\nabla}}{\mu}\cdot\boldsymbol{\mathbf{n}}=\dot{m}(\boldsymbol{\mathbf{x}})$, where $\dot{m}(\boldsymbol{\mathbf{x}})$ is the desired mass flux. In particular, we will set $\dot{m}(\boldsymbol{\mathbf{x}})=0$ at the solid boundary, $\Gamma_{S}$. In summary, the equations we are solving in this study are the following:

[TABLE]

In addition to the conservation of mass, Eqs. 5a, 5b, 5c and 5d imply the following energy-conservation law, involving the phase field and the chemical potential:

[TABLE]

An advantage of the cubic surface energy (1) over other surface energy formulations (see [33] for a review of wetting boundary conditions for binary fluids) is that the well-known hyperbolic tangent profile is an equilibrium solution in more than 1 dimension. Specifically, the function

[TABLE]

is solution to the CH equation posed in the half plane $\{y\geq 0\}$ with the boundary condition (1) at $\{y=0\}$ and constant $\theta(\boldsymbol{\mathbf{x}})=\theta$. A schematic representation of this solution and the corresponding fluid-fluid interface is given in Fig. 1.

A drawback of the cubic wall energy (1) is that the conservation of energy no longer seems to imply stability bounds for the solution, making it impossible to use the tools traditionally employed (see e.g. [20]) to prove the existence of a solution. Indeed, an application of the trace inequality gives only that, under appropriate regularity assumptions on $\phi$:

[TABLE]

where we used Hölder’s inequality and Young’s inequality with a parameter. Therefore, the wall energy cannot be controlled by the mixing energy for arbitrary domains. This issue can be remedied by a simple modification of the wall energy outside of the physical range $[-1;1]$; instead of (1), we consider the following wall energy:

[TABLE]

This function is such that $F_{w}^{*}(\phi)=F_{w}(\phi)$ for $\phi\in[-1,1]$, $F_{w}^{*}\in C^{2}()$, and $(F_{w}^{*})^{\prime\prime}$ is absolutely continuous, which makes it possible to prove the second order convergence of our time-stepping scheme, see Section 3. Another possibility would have been to choose constant values for $F_{w}^{*}$ outside of the interval $[-1,1]$, but this would have lead to $F_{w}^{*}$ being only $C^{1}()$, making it more difficult to show second order convergence theoretically. The weak formulation of Eqs. 5a, 5b, 5c and 5d with the modified wall energy (9) is as follows: find $(\phi,\mu)$ such that

[TABLE]

and the following variational formulation is satisfied:

[TABLE]

with $f_{w}^{*}:=(F_{w}^{*})^{\prime}$ and where $\langle\cdot,\cdot\rangle$, $(\cdot,\cdot)$ and $(\cdot,\cdot)_{\partial\Omega}$ denote respectively the duality pairing between $(H^{1}(\Omega))^{\prime}$ and $H^{1}(\Omega)$, the standard inner product in $L^{2}({\Omega})$, and the standard inner product in $L^{2}({\partial\Omega})$. For simplicity of notations, the symbols $F_{w}$, $f_{w}$ and $E$ will refer in the rest of this paper to $F_{w}^{*}$, $f_{w}^{*}$, and $E_{m}+\int_{\partial\Omega}F_{w}^{*}\,\mathrm{d}\sigma$, respectively.

2.1 Existence of a solution

We can show the following existence result for the weak formulation of the Cahn–Hilliard system with the modified boundary condition presented above, under appropriate regularity assumptions for the initial condition and the mass flux $\dot{m}$.

Theorem 2.1**.**

Assume that $\phi_{0}\in H^{1}(\Omega)$ and $\dot{m}\in C([0,T];L^{2}({\partial\Omega}))$. Then there exists a pair of functions $(\phi,\mu)$, with

$\phi\in L^{\infty}(0,T;H^{1}(\Omega))\bigcap C([0,T];L^{2}({\Omega}))$, 2. 2.

$\mathchoice{\frac{\partial{\phi}}{\partial t}}{\partial_{t}\phi}{}{}\in L^{2}({0,T;(H^{1}(\Omega))^{\prime}})$, 3. 3.

$\phi(0)=\phi_{0}$, 4. 4.

$\mu\in L^{2}({0,T;H^{1}(\Omega)})$,

*that solve the variational formulation Eqs. 10a and 10b. *

Proof.

See Appendix A. ∎

3 Numerical method

In this section we introduce a new time-stepping scheme to solve the CH equation (3) with the non-linear wetting boundary condition (2), which is a generalization of the optimal dissipation scheme of order 2, OD2, developed in [27]. We decided to extend this particular scheme because the authors of [27] showed that, among all the linear schemes they proposed, it is the most accurate and the least dissipative. And in selected test cases, they showed that for a large enough time step, it is the only scheme that leads to the correct equilibrium solution. We refer to our scheme as OD2-W, with W denoting wetting, and show that it leads to a consistent discrete energy law.

We also develop a new adaptive time-stepping strategy which, combined with adaptation in space, leads to a fully adaptive finite element method. An excellent introduction to the finite element method and corresponding mixed formulations can be found in [13] and to mesh generation and adaptive refinement in [24].

3.1 OD2-W scheme

In this section, we assume for simplicity that $\dot{m}=0$ and that $\theta$ is uniform on $\partial\Omega$. We denote by $\Delta t$ the time step, and by $\phi^{n}$ and $\mu^{n+1/\alpha}$ the numerical approximations of $\phi$ and $\mu$ at times $t^{n}$ and $t^{n}+\frac{1}{\alpha}\Delta t$, respectively. To define a discretization in time of the CH system appropriate for wetting phenomena, we follow the approach proposed in [27] to design an optimal dissipation scheme, and consider the following generic implicit-explicit numerical scheme: given $\phi^{n}\in H^{1}(\Omega)$, find $(\phi^{n+1},\mu^{n+\frac{1}{\alpha}})\in H^{1}(\Omega)\times H^{1}(\Omega)$ such that,

[TABLE]

In these expressions, $\hat{f}_{m},\hat{f}_{w}$ are functions to be specified, linear in their second argument. The parameter $\alpha\in\{1,2\}$ determines the accuracy of the numerical scheme, and the parameter $\beta\in[0,1-1/\alpha]$ controls the numerical diffusion. The function $\phi^{n+\frac{1}{\alpha}+\beta}$ is defined by linear interpolation between $\phi^{n}$ and $\phi^{n+1}$,

[TABLE]

and $\delta_{t}\phi^{n+1}$ is the approximation of the time derivative of $\phi$ given by

[TABLE]

In most numerical experiments presented in this chapter, we consider the case $(\alpha,\beta)=(2,0)$ (OD2-W), but we note that other usual choices include $(\alpha,\beta)=(1,0)$ (OD1-W) and $(\alpha,\beta)=(2,\mathcal{O}(\Delta t))$ (OD2mod-W). By taking $\psi=\mu^{n+\frac{1}{\alpha}}$ and $\nu=\delta_{t}\phi^{n+1}$ in (11), we obtain

[TABLE]

where $ND(\phi^{n},\phi^{n+1})$, representing the non-physical numerical dissipation introduced by the time-stepping scheme, can be broken down in three parts:

[TABLE]

with

[TABLE]

Notice that the philic dissipation is always nonnegative, with $ND_{\textrm{philic}}(\cdot,\cdot)=0$ if $(\alpha,\beta)=(2,0)$ (OD2-W), $ND_{\textrm{philic}}(\cdot,\cdot)=\mathcal{O}(\Delta t^{2})$ if $(\alpha,\beta)=(2,\mathcal{O}(\Delta t))$ (OD2mod-W), and $ND_{\textrm{philic}}(\cdot,\cdot)=\mathcal{O}(\Delta t)$ if $(\alpha,\beta)=(1,0)$ (OD1-W). The two other terms can be expanded using Taylor’s formula, taking into account that $F_{m}$ is a polynomial of degree 4 and using the integral form of the remainder:

[TABLE]

This suggests the following choices for the functions $\hat{f}_{m}$ and $\hat{f}_{w}$:

[TABLE]

where the last expression is convenient for programming purposes. We note that this methodology to derive a second-order scheme can be applied mutatis mutandis when using the unmodified wall energy (1), although we haven’t been able to prove that the weak formulation admits a solution in that case. Doing so leads to $\hat{f}_{w}(\phi^{n},\phi^{n+1})=-(\sqrt{2}/2)\,\cos{\theta}\,(1-\phi^{n}\,\phi^{n+1})$, which coincides with (15b) when $\phi^{n}\in[-1,1]$. In either case, we have the following property:

Property 3.1**.**

Assume that $\alpha=2$ and $\beta=0$. Then the numerical dissipation term in Eq. 12 is such that

[TABLE]

with $C:=\left(C_{1}\,\|\phi\|_{C([0,T],L^{\infty}(\Omega))}\,\|\partial_{t}\phi\|_{C([0,T],L^{3}(\Omega))}^{3}+C_{2}\,\|\partial_{t}\phi\,\|_{C([0,T],L^{3}(\partial\Omega))}^{3}\right)$, provided that all the terms in the definition of $C$ are well-defined.

Proof.

In [27], the authors show that:

[TABLE]

For the wall term, we obtain from Eqs. 14b and 15b:

[TABLE]

∎

In addition to the energy law (12), the numerical scheme (11) satisfies a discrete version of the conservation law (4) presented in Section 2, which can be seen by choosing $\psi=1$ in Eq. 11a.

Property 3.2**.**

The numerical solution satisfies the following mass conservation law:

[TABLE]

3.1.1 Space discretization and adaptive mesh refinement

Our approach for mesh adaptation is based on a method proposed in [30], and implemented through the FreeFem++ functions adaptmesh (in 2D) and mshmet (in 3D). The idea of the method is to define a metric on the computational domain based on the solution at the current time step, and to use for the next time step a mesh that is uniform in that metric. The metric we consider corresponds to the following metric tensor, depending only on the phase field $\phi$:

[TABLE]

where $(\lambda_{i}(\boldsymbol{\mathbf{x}}))_{i=1}^{d}$ are the eigenvalues of the Hessian of $\phi$ at $\boldsymbol{\mathbf{x}}$, $R(\boldsymbol{\mathbf{x}})$ is the matrix containing the associated orthonormal eigenvectors, and $\gamma>0$ is a parameter controlling the interpolation error. A standard algorithm of Delaunay type is used to generate a mesh that is equilateral and uniform with characteristic length 1 in that metric. This mesh definition ensures that the interpolation error of the phase field is roughly equi-distributed over the parts of the domain where $h_{\max}^{-2}\leq\frac{1}{\gamma}\max_{i\in\{1,2,\dotsc,d\}}\left|\lambda_{i}\right|\leq h_{\min}^{-2}$.

In most of the simulations presented in the next section, we set $h_{\min}$ to a value lower than or equal to $\varepsilon/5$, to ensure that enough mesh points are available for the discretization of the interface region in its normal direction, and $h_{\max}$ to a value small enough that a good approximation of the chemical potential is possible. For 3D simulations, however, choosing $h_{\min}\leq\varepsilon/5$ when $\varepsilon$ is of the order of $0.01$ leads to a prohibitive computational cost; in these cases we have thus used a less precise value, as specified in the relevant sections.

For a given mesh \mbox{\cal{T}}=\bigcup_{i=1}^{N_{T}}T_{i}, we use the standard finite element space

[TABLE]

with $P_{\rho}$ the space of polynomials of degree $\rho$. In the numerical experiments below, we used both quadratic elements ($\rho=2$) and linear ones ($\rho=1$). Space discretization is achieved by replacing $H^{1}(\Omega)$ by $V_{h}$ in the variational formulation (11), leading to a sparse unsymmetric linear system at each iteration.

3.1.2 Time step adaptation

Here we assume that $\dot{m}=0$ in the boundary condition (5d). From Eqs. 4 and 6, this implies that $M(\phi)$ is constant in time and $E(\phi)$ decreases. Numerical exploration suggests that large free-energy variations are usually caused by topological changes of interfaces, corresponding to physical phenomena such as the coalescence of droplets. Since capturing such phenomena precisely is crucial to the accuracy of the solution, we propose an adaptive strategy that aims to limit the variation of free energy at each time step. We adapt the time step based on the dissipation of free energy,

[TABLE]

which is equal to $E(\phi^{n+1})-E(\phi^{n})$ up to numerical dissipation. Here $\Delta t^{n+1}:=t^{n+1}-t^{n}$. Six parameters enter in our time-adaptation scheme:

• •

$\Delta t_{\min},\Delta t_{\max}$: the time steps below which we stop refining and beyond which we stop coarsening, respectively.

• •

$\Delta E_{\min}$: the variation of free energy below which we increase the time step at the next iteration.

• •

$\Delta E_{\max}$: the variation of free energy beyond which we refine the time step and recalculate the numerical solution.

• •

$f>1$: the factor by which the time step is multiplied or divided at each adaptation.

• •

$0<\zeta\ll 1$: a small parameter controlling the maximum amount by which the free energy is allowed to increase at each time step. In the numerical experiments presented in Section 4, we always take $\zeta=0.01$.

The condition $(E(\phi^{*})-E(\phi^{n}))>\zeta\,\Delta E_{\max}$ serves to guarantee that the method does not blow up. The choice of a nonzero right-hand side is motivated by the fact that, when the system is close to equilibrium, it can happen that $E(\phi^{*})>E(\phi^{n})$. This is because, in contrast with the sign of $ND_{\textrm{philic}}(\phi^{n},\phi^{n+1})$, which is always positive or zero according to Eq. 13, the signs of $ND_{\textrm{phobic}}(\phi^{n},\phi^{n+1})$ and $ND_{\textrm{wall}}(\phi^{n},\phi^{n+1})$ are in general unknown.

In the numerical experiments presented in Section 4, we chose $\Delta t_{\min}=0$. Since the numerical dissipation term scales as $\Delta t^{2}$ for the method OD2-W, the inequality $E(\phi^{n+1})\leq E(\phi^{n})+\zeta\,\Delta E_{\max}$ will always hold for $\Delta t$ small enough, so the refinement process is guaranteed to terminate at each iteration.

If the condition $\dot{m}=0$ is not satisfied, then the adaptation scheme we proposed, being based on limiting the decrease and preventing the increase (beyond a very small tolerance equal to $\zeta\Delta E_{\textrm{max}}$) of free energy at each time step, is not applicable. Indeed, by Eq. 6, a nonzero mass flux can cause a positive variation of the total free energy, which our scheme would be unable to distinguish from a nonphysical increase of energy caused by numerical errors. In this nonzero flux situation, the adaptation scheme presented in [28], which is based on limiting the numerical dissipation, could be used instead (with suitable modifications to accommodate for the presence of a flux).

4 Numerical results

The new numerical method is applied on a number of test cases. We have used FreeFem++ [29] for the implementation off the finite element method and 2D mesh adaptation, umfpack [16] for the linear solver, mshmet [23] and tetgen [53] for the mesh adaptation in 3D, and gmsh [26] for the description of the geometry, post-processing and 3D visualization. In Section 4.1 we check that the numerical scheme leads to the correct equilibrium solution in the simple case of a droplet spreading on a philic or phobic substrate. In Section 4.2 we study the convergence of the method with respect to the time step and the mesh size, when a uniform mesh and a constant time step are used. In Section 4.2 we illustrate the time-adaptation scheme in the case of two droplets coalescing on a substrate. Finally, Section 4.4 demonstrates the ability of the numerical scheme to scrutinize wetting phenomena in more complicated geometries, and in the presence of heterogeneous substrates. The code used for the simulations is available on line, see [62].

4.1 Equilibrium contact angle

We consider a 2D sessile droplet on a flat substrate where we impose the no-flux condition and the wetting condition (2) with the modified wall energy (9) and uniform contact angle $\theta$:

[TABLE]

Our aim in this section is to check that our method is able to accurately capture the imposed contact angle, $\theta$. Figure 2 shows the equilibrium position of a droplet for different values of $\theta$, for $b=1$ and $\varepsilon=5\times 10^{-3}$. In all cases we used the scheme OD2-W, with quadratic basis functions and adaptation in space only using the parameters $h_{\max}=10\,h_{\min}=0.01$, and we computed the contact angle of the $\phi=0$ isoline at the substrate. A very good agreement is achieved between the imposed equilibrium contact angle and the observed numerical one.

4.2 Convergence of the method

Here, we study the convergence of the method when both time step and mesh size decrease. The problem we considered to that purpose is the coalescence of two adjacent sessile droplets as they spread on a flat substrate. For the simulation, we used the initial condition

[TABLE]

in the domain $[0,2]\times[0,0.5]$, with $x_{1}=0.65$, $x_{2}=1.35$, $r=0.25$, and at the boundary we imposed a uniform contact angle, $\theta=\pi/4$, using the wall energy (9). Only linear elements were used.

For the convergence as $h\to 0$, we solved the problem numerically for several values of $h$, without mesh adaptation and for $\varepsilon=0.1$, so that enough data points could be generated at a reasonable numerical cost. Since the exact solution to the CH equation in this case is not known analytically, we calculated the error by comparison of the numerical solutions to the solution obtained with the smallest value of $h$. Results are presented in Fig. 3. As we can see, the observed convergence rate is almost equal to 2, which is the optimal rate in the case of linear basis functions.

Now we address the convergence with respect to the time step. For this case, we used the parameters $\varepsilon=0.1$, $b=10^{4}$, and the minimum time step we considered was $\Delta t^{*}:=0.00665.$ In Fig. 4, we present convergence curves for OD1-W, OD2-W, and OD2mod-W. We note that the convergence rates are close to the expected ones, and that the use of OD2 gives significantly more accurate results than the other two methods. In Fig. 5, the total numerical dissipation produced by the numerical schemes is presented. Here too, numerical results agree with the theoretical results of Section 3.

4.3 Time-adaptation scheme

In this section, we examine the performance of the adaptive time-stepping scheme in the case of two droplets evolving on a chemically homogeneous substrate. We start from the situation where $\phi=-1$ everywhere except in two half-circles, of radius $r=0.25$ and centred at $(0.65,0)$ and $(1.35,0)$, where $\phi=1$. We used quadratic basis functions and the following parameters: $b=10^{-4}$, $\varepsilon=0.01$, $f=\sqrt{2}$, $\Delta t_{0}=0.02$, $\Delta t_{\min}=0$, $\Delta t_{\max}=16\Delta t_{0}$, $\Delta E_{\min}=0.0001$, $\Delta E_{\max}=0.0002$, $h_{\max}=0.05$ $h_{\min}=0.001$, and for $\theta$ we considered three values: $\pi/4$, $\pi/2$, $3\pi/4$.

Snapshots of the phase field and of the chemical potential at different times of the simulation are presented in Figs. 6 and 7 for the case $\theta=\pi/4$ and $\theta=3\pi/4$, respectively. The case $\theta=\pi/2$ is less interesting because, in view of the initial condition, the droplets remain essentially motionless throughout the simulation; we do not present snapshots of the solution in that case.

The evolution of the time step, of the number of recalculations, and of the free energies is displayed in Fig. 8. In all three cases, the time step is refined several times at the first iteration, which can be explained as follows. On one hand, this refinement originates from the discontinuity of the initial condition at the interface between the two fluids. On the other hand, in the cases $\theta=\pi/4$ and $\theta=3\pi/4$, it is also a consequence of the fact that the initial condition does not satisfy the wetting boundary condition. Since the initial angle between the interface and the substrate is equal $\pi/2$, the number of recalculations performed at the first iteration is higher for $\theta=\pi/4,3\pi/4$ than for $\theta=\pi/2$. After the initial refinement, the time step steadily increases to its maximum allowed value for $\theta=\pi/2$ and $\theta=3\pi/4$, but when $\theta=\pi/4$ a second refinement occurs to capture the coalescence of the droplets.

In this latter case, we observe, simultaneous with the second refinement of the time step, an increase in the rate of dissipation of free energy. After the formation of a new stable interface, the total free energy continues to decrease, but more slowly, as a new droplet, formed by the merging of the two original droplets, moves towards its equilibrium position. We clearly identify the coalescence time by looking at the singularity in the curve corresponding to the mixing energy. This energy increases before coalescence, as the interfaces are being stretched, and it decreases steadily after. The wall energy, on the other hand, decreases at first and increases in the later stage of the simulation. As prescribed by Algorithm 1, the time step detects the variations of free energy; it decreases when the rate of variation of the total free energy increases, and conversely.

For comparison purposes, we also included in Figs. 8, 8 and 8 data corresponding to the case where a fixed time step is used for the simulations presented in this section. There does not currently exist any result with conditions on the time step that ensure the stability of OD2, and we haven’t been able to show stability results for OD2-W either. In practice, we observed that the time step required to ensure stability of OD2-W with the set of parameters we use in this test case would lead to a very high computational cost. We point out that, contrary to what we expected, the time step required to achieve stable integration in time with the modified wall energy (9), which we use in this paper, seems to be generally smaller than with the cubic formulation (1). To keep the computational cost at a reasonable level, we carried out the simulations with a fixed time step using the method OD1-W, the greater stability of which enabled us to choose $\Delta t=0.02$. In Figs. 8, 8 and 8, we see that, for the same contact angle, the curves corresponding to a fixed and an adaptive time step are almost undistinguishable. The agreement is also very good at the level of the phase field and chemical potential, although we do not present snapshots of the solutions obtained with a fixed time step.

The CPU times corresponding to the three contact angles considered are presented in Table 1. Without adaptation, the simulations take significantly longer to run, which is consistent with the fact that more iterations (20000) were necessary to reach the final time. In addition, among the simulations that used an adaptive time-step, the difference between the CPU times is also significant, with the case $\theta=\pi/4$ taking more than twice as long as the case $\theta=\pi/2$.

4.4 Wetting in complex geometries and with heterogeneous substrates

We now present the results of numerical experiments in more complicated and realistic settings, in both 2D and 3D systems. The results presented in this section were all obtained using linear basis functions.

4.4.1 3D droplet on a chemically heterogeneous substrate

We study the dynamics of a 3D sessile droplet on a flat substrate with chemical heterogeneities, i.e. the contact angle has a spatial dependence now, say $\theta=\theta(x,y)$. This situation typically arises in electro-wetting settings [42]. It is widely accepted that the droplet shape can be controlled using patterned substrates, e.g. Ref [63, 52], that may also be modelled efficiently using a space varying contact angle [63]. We consider chemical heterogeneities on the substrate of the form

[TABLE]

with $\theta_{0}=\frac{\pi}{2}$ the mean contact angle, $a=\frac{\pi}{6}$ the amplitude, and $f_{x}=f_{y}=4$ the frequencies in $x$ and $y$ directions, respectively. As initial condition we take a droplet of base radius $r_{0}=0.24$ centered at $\boldsymbol{\mathbf{x}_{0}}=(0.5,0.5,0)$. The initial values of the phase field are given as

[TABLE]

Results are displayed in Fig. 9. The droplet, initially spherical, spreads on the hydrophilic regions of the substrate, and retracts from the hydrophobic patches. While we do not present any quantitative analysis of the error in this case, we note that the wetting behaviour agrees qualitatively with what one might expect intuitively from our understanding of wetting phenomena. While it progresses towards equilibrium, the droplet adopts a diamond-like shape.

For this test case, we used the method OD2-W with adaptation in space and time. The parameters used were the following: $b=10^{4}$, $\varepsilon=0.02$, $h_{\max}=10\,h_{\min}=0.1$, $\Delta t_{0}=0.0016$, $\Delta t_{\min}=0$, $\Delta t_{\max}=16\,\Delta t_{0}$, $f=\sqrt{2}$, $\Delta E_{\max}=10\,\Delta E_{\min}=0.0001$. With these parameters, the time step was refined only at the beginning of the simulation, which is consistent with the absence of coalescence events in this case. There were 24 recalculations at the first time step, corresponding to a refinement of the time step by a factor $f^{24}=4096$.

4.4.2 Diffusion in a 3D porous medium

Here we consider a binary fluid in a model porous medium consisting of a cube filled with spheres. The cube has edges of length 1, and the spheres have radius $0.1$ and are located at positions $(1.5,1.5,1.5)+2\Delta(i,j,k)$ with $\Delta=1/7$ and $i,j,k\in\{0,1,2\}$. We take all the substrates to be neutral, i.e. $\theta=\frac{\pi}{2}$, and the initial condition is the same as used before, defined by Eq. 19. In addition, we include an inflow boundary condition at the bottom of the cube to represent a pore where liquid can be pumped in. The radius of this pore is $0.1$ and is located at $(0.5,0.5,0)$. This boundary condition can be incorporated by imposing

[TABLE]

which models the situation when the component $\phi=1$ is pumped into the domain. Under these conditions, we study how the flow is affected by the geometry of the domain. Our results are depicted in Figs. 10 and 11.

The imposed contact angle at the spheres is $\pi/2$, forcing the isosurface to stay normal to the spheres as long as these are not completely covered. Because of the boundary condition (20), the mass increases linearly, and the free energy increases, in agreement with Eqs. 4 and 6. This case study demonstrates the ability of our method to easily tackle complex geometries. The parameters used for this test case are the same as in Section 4.4.1, except that we employed the fixed time step $\Delta t=0.001$.

4.4.3 Nucleation processes with complex boundaries

The last problem we study is the process of phase separation in a domain with complex boundary characterized by different length scales. Specifically, we consider a domain defined by the coastline of the two islands that form the United Kingdom and Ireland. Starting from a satellite black and white picture, we extracted the isolines that define the contour of the different islands, which we passed to the FreeFem++ mesh generator to obtain a triangular mesh (for this, we based our code on a FreeFem++ example for the Leman lake). At the boundary we consider the contact angles $\theta=\pi/4,\pi/2,3\pi/4$, and we assume that the phase field is initially set to a random value at each grid point, drawn from a random normal distribution with variance $0.1$. A fixed mesh was used for this simulation, and the parameters used were $b=1000$, $\varepsilon=0.02$, $\Delta E_{\min}=0.02$, $\Delta E_{\max}=0.04$, $f=\sqrt{2}$, $\Delta t_{\min}=0$, $\Delta t_{\max}=1$.

The evolution of the phase field and of the chemical potential in the case $\theta=\pi/4$, obtained with the adaptive time-stepping scheme 1, is presented in Fig. 12. For each of the contact angles considered, we also ran a simulation with the fixed time step $\Delta t=0.01$, using the method OD1-W instead of OD2-W to benefit from the stabilizing effect introduced by the philic numerical dissipation of OD1-W. We note in particular that OD2-W is unstable for the selected value of $\Delta t$, with oscillations appearing in the energy curves from the first iterations, and that the time step would have to be reduced significantly to ensure stability. The final configurations (time 500) are presented in Fig. 13 for the three contact angles considered. We observe that the final configurations are different depending on whether or not an adaptive time step is used, which can be attributed to the high sensitivity of the solution to perturbations of the initial condition chosen for this test case; the areas where separation of the phases first occurs is influenced by numerical errors in the early stages of the simulation.

Simulation data are presented in Fig. 14. With an adaptive time-stepping scheme, it appears from Fig. 14 (a) that, overall, the time step increases steadily as the frequency of coalescence events decreases. At specific times, the time step decreases slightly in order to accurately capture the evolution. As expected, the total free energy has a roughly constant negative slope when plotted against the iteration number. Here too, we observe a small discrepancy between the fixed and adaptive cases, which is consistent with differences observed at the final time in Fig. 13.

The CPU times corresponding to the simulations presented in this section are displayed in 2. For the parameters selected, the adaptive time-stepping scheme leads to a lower computational cost.

This test demonstrates the advantage of using a finite element approach, as it would have been very complicated to solve the CH equation in the geometries we consider here with e.g. a spectral method or finite differences.

5 Conclusions

We have proposed a new, fast and reliable numerical method to solve the CH equation with a wetting boundary condition. Our method is a generalization of the OD2 scheme introduced in [27], which considered only the homogeneous condition $\boldsymbol{\mathbf{\nabla}}\phi\cdot\boldsymbol{\mathbf{n}}=0$. In addition, we have designed a new time-step adaptation algorithm, leading to a scheme that is adaptive both in space and time, and we have shown that this scheme is mass-conservative and satisfies a consistent discrete energy law.

We checked the validity of the proposed numerical scheme with several examples. First we considered the relaxation towards equilibrium of a sessile droplet and the coalescence of two sessile droplets on flat, chemically homogeneous substrates; then we considered several multiphase systems in complex geometries or surrounded by chemically heterogeneous substrates.

Compared to finite differences or spectral approaches, the method introduced here has the advantage that it can be used without modification with complex geometries. Furthermore, the numerical scheme we have proposed can easily be extended to include at least two additional features. First, a linear, energy-stable, second-order scheme could be developed for the three-component CH model with wetting boundary conditions, building on the work of [6, 7].

Second, we remark that in our work, we considered a regime in which contact line motion is controlled by diffusive interfacial fluxes, or in other words, we considered a large diffusivity limit, where any possible advection effects are neglected. To account for such effects the model must be appropriately modified to include an advection term coupled to the Navier–Stokes equations [1, 51, 35, 36, 37, 65]. Such generalizations are indeed possible within the proposed numerical scheme and we hope to address these and related issues in future studies.

6 Acknowledgements

We are grateful to Ben Goddard, Demetrios Papageorgiou and Srikanth Ravipati for useful suggestions. We acknowledge financial support by the UK Engineering and Physical Sciences Research Council (EPSRC) through Grants No. EP/L027186 and EP/L020564 and European Research Council (ERC) via Advanced Grant No. 247031.

Appendix A Proof of Theorem 2.1

Before presenting the proof, we recall a particular Sobolev embedding for smooth bounded domains; see e.g. [22, Chap. 5]. Let $d\geq 2$, $\emptyset\neq\Omega\subset^{d}$ be open with $C^{1}$ boundary, and assume that $q<\infty$ if $d=2$ or $q<p^{*}:=2d/(d-2)$ if $d>2$. Then the following embedding is compact

[TABLE]

We also recall two other well-known compactness results; see e.g. [41]. Let $X,Y,Z$ be Banach spaces with a compact embedding $X\subset Y$ and a continuous embedding $Y\subset Z$. Then the following embeddings are compact:

[TABLE]

Proof.

Without loss of generality, we assume that the mobility, $b$, is equal to $1$. In the spirit of [20, Theorem 2], we apply a Faedo–Galerkin approximation. Let $\{\varphi_{n}\}_{n\in\mathbb{N}}$ and $\{\lambda_{n}\}_{n\in\mathbb{N}}$ denote the eigenfunctions and eigenvalues of the Laplace operator with a homogeneous Neumann boundary condition, i.e.

[TABLE]

normalized such that

[TABLE]

We assume without loss of generality that $\lambda_{1}=0$. To build an approximation of the solution to Eqs. 10a and 10b in the finite-dimensional space $S_{N}:=\text{span}\{\{\}\varphi_{1},\dots,\varphi_{N}\}$, we consider the following ansatz,

[TABLE]

and the variational formulation

[TABLE]

To this formulation corresponds the following system of ordinary differential equations, with unknown functions ${\left\{a^{N}_{n}\right\}}_{n=1}^{N}$ and ${\left\{b^{N}_{n}\right\}}_{n=1}^{N}$:

[TABLE]

for $n=1,\dots,N$. Local existence and uniqueness of a solution to this system of equations is guaranteed by the fact that that the right-hand side of (23a) depends continuously on the coefficients ${\left\{a^{N}_{n}\right\}}_{n=1}^{N}$. To show the existence of a global solution, we will use the a priori estimate presented in the following lemma.

Lemma A.1**.**

Assume that $\left|F_{w}(\phi)\right|\leq C(1+\left|\phi\right|^{2})$. Then the solution $(\phi^{N},\mu^{N})$ to Eqs. 22a, 22b and 22c satisfies

[TABLE]

where $C$ is independent of $N$ and $\Omega_{T}:=\Omega\times(0,T)$.

Proof.

Setting $\bar{\phi}=\mu^{N}$, $\bar{\mu}=\partial_{t}\phi^{N}$ in  Eqs. 22a and 22b and subtracting leads to the equation

[TABLE]

Using a trace inequality, Hölder’s inequality, and Young’s inequality with a parameter, we have, for all $u\in H^{1}(\Omega)$,

[TABLE]

Now we use the simple fact that, for any $\beta>0$ and $0\leq s\leq t$, the inequality $\left|x\right|^{s}\leq\beta^{s}+\beta^{s-t}\,\left|x\right|^{t}$ holds true for all $x\in$, to obtain

[TABLE]

for a constant $C$ independent of $u$.

In addition, using a trace inequality, Poincaré inequality, and  (22b) with $\bar{\mu}=1$,

[TABLE]

Integrating (25) in time, and rearranging using Eqs. 26 and 27,

[TABLE]

where we used the notations $\Omega_{t}$ and $\partial\Omega_{t}$, $t>0$, to denote $\Omega\times(0,t)$ and $\partial\Omega\times(0,t)$, respectively. The last inequality holds by the assumptions that $\phi_{0}\in H^{1}(\Omega)$ and $\dot{m}\in C([0,T];L^{2}({\partial\Omega}))$. Using a Grönwall inequality, we have Eq. 24. ∎

By integration by parts of the first term in Eq. 24, we obtain $\sum_{n=1}^{N}\lambda_{n}(a^{N}_{n})^{2}<C$. This result, together with the inequality

[TABLE]

implied by Equation 23a and the fact that $\lambda_{1}=0$, show that the coefficients ${\left\{a^{N}_{n}\right\}}_{n=1}^{N}$ do not blow up, and by Eq. 23b neither do the coefficients ${\left\{b^{N}_{n}\right\}}_{n=1}^{N}$, implying global existence.

In addition to (24), we have the usual estimate on $\mathchoice{\frac{\partial{\phi^{N}}}{\partial t}}{\partial_{t}\phi^{N}}{}{}$: denoting by $\Pi^{N}$ the $L^{2}(\Omega)$ projection on $S_{N}$, for all $\psi\in L^{2}({0,T;H^{1}(\Omega)})$ the following holds:

[TABLE]

This shows that $\left\|\mathchoice{\frac{\partial{\phi^{N}}}{\partial t}}{\partial_{t}\phi^{N}}{}{}\right\|_{L^{2}({0,T;(H^{1}(\Omega))^{\prime}})}\leq C$.

Let $p$ be such that the embedding $H^{1}(\Omega)\subset L^{p}({\Omega})$ is compact, i.e., by Rellich–Kondrachov theorem, $p<\infty$ if $d=1$ or $d=2$, and $p<\frac{2d}{d-2}$ if $d>2$. Using Eqs. 24 and 28, we can apply results (21a) and (21b) to our case, with $X=H^{1}(\Omega)$, $Y=L^{p}({\Omega})$ and $Z=(H^{1}(\Omega))^{\prime}$, to conclude that there exists a subsequence such that

[TABLE]

when $N\to\infty$. In addition, note that since $\phi^{N}$ is bounded in $L^{\infty}({0,T;L^{2}({\partial\Omega})})$, there is a subsequence such that $\phi^{N}\to v$ weak-* in $L^{\infty}({0,T;L^{2}({\partial\Omega})})$ for some function $v$ in that space, and thus also $\phi^{N}\to v$ weakly in the coarser $L^{2}({0,T;L^{2}({\partial\Omega})})$. But also $\phi^{N}\to\phi$ weakly in $L^{2}({0,T;L^{2}({\partial\Omega})})$, because $\phi^{N}\to\phi$ weakly in $L^{2}({0,T;H^{1}(\Omega)})$ and by continuity of the trace operator (indeed, an operator between Hilbert spaces that is continuous in the strong topologies, is continuous in the weak ones too), so $v=\phi$. The same reasoning can be applied to $f_{w}(\phi^{N})$, taking into account that $f_{w}:u\rightarrow f_{w}(u)$ is continuous on $L^{2}({0,T;L^{2}({\partial\Omega})})$, to conclude

[TABLE]

Regarding the chemical potential, testing (22b) with $\bar{\mu}=1$ implies that:

[TABLE]

which, together with the energy estimate (24), implies that $\mu^{N}$ is bounded in $L^{2}({0,T;H^{1}(\Omega)})$, leading to the existence of a further subsequence such that

[TABLE]

Proceeding in a standard fashion, we consider an integer $M$ and arbitrary functions $\phi^{M},\mu^{M}\in C([0,T],H^{1}(\Omega))$ such that

[TABLE]

with $\{\bar{a}_{n}\}_{n=1}^{M}$, $\{\bar{b}_{n}\}_{n=1}^{M}$ smooth functions. Using $\phi^{M}$ and $\mu^{M}$ as test functions in Eqs. 22a and 22b, integrating in time, taking the limit $N\to\infty$, and using the convergence results given in Eqs. 29b, 29c, 29d, 30, 31 and 29a, we obtain

[TABLE]

from which we conclude using a standard density argument. ∎