An Integral Equation Method for the Cahn-Hilliard Equation in the Wetting Problem

TL;DR

This paper introduces an efficient integral equation method for solving the Cahn-Hilliard equation with boundary conditions modeling wetting, utilizing advanced numerical techniques for high accuracy and linear complexity.

Contribution

It develops a novel integral equation approach combining volume potentials, SKIE, and QBX with FMM acceleration for the Cahn-Hilliard equation in wetting problems.

Findings

Achieves high-order convergence with adaptive refinement.

Linear complexity in degrees of freedom.

Effective handling of boundary conditions with prescribed Young's angles.

Abstract

We present an integral equation approach to solving the Cahn-Hilliard equation equipped with boundary conditions that model solid surfaces with prescribed Young's angles. The discretization of the system in time using convex splitting leads to a modified biharmonic equation at each time step. To solve it, we split the solution into a volume potential computed with free space kernels, plus the solution to a second kind integral equation (SKIE). The volume potential is evaluated with the help of a box-based volume-FMM method. For non-box domains, source density is extended by solving a biharmonic Dirichlet problem. The near-singular boundary integrals are computed using quadrature by expansion (QBX) with FMM acceleration. Our method has linear complexity in the number of surface/volume degrees of freedom and can achieve high order convergence with adaptive refinement to manage error from function extension.

Figures

Tables

Table 1. Table 1: Results Using First-Order Approximation ( q b subscript 𝑞 𝑏 q_{b} : number of quadrature points per boundary cell, q v subscript 𝑞 𝑣 q_{v} : number of quadrature points per volume cell, e 0 subscript 𝑒 0 e_{0} : relative error at t = δ ​ t 𝑡 𝛿 𝑡 t=\delta t ).

$q_b$	$q_v$	$\delta x$	$e_0$	EOC
1	1	$1.56\times {10}^{-2}$	$9.91\times {10}^{-1}$	–
1	1	$7.81\times {10}^{-3}$	$6.72\times {10}^{-1}$	0.56
1	1	$3.91\times {10}^{-3}$	$2.91\times {10}^{-1}$	1.21
1	1	$1.95\times {10}^{-3}$	$1.34\times {10}^{-1}$	1.19
1	1	$9.77\times {10}^{-4}$	$9.02\times {10}^{-2}$	0.57
1	1	$4.88\times {10}^{-4}$	$5.24\times {10}^{-2}$	0.78

Table 2. Table 2: Results Using Second-Order Approximation ( q b subscript 𝑞 𝑏 q_{b} : number of quadrature points per boundary cell, q v subscript 𝑞 𝑣 q_{v} : number of quadrature points per volume cell, e 0 subscript 𝑒 0 e_{0} : relative error at t = δ ​ t 𝑡 𝛿 𝑡 t=\delta t ).

$q_b$	$q_v$	$\delta x$	$e_0$	EOC
2	4	$1.56\times {10}^{-2}$	$4.01\times {10}^{-1}$	–
2	4	$7.81\times {10}^{-3}$	$2.41\times {10}^{-1}$	0.73
2	4	$3.91\times {10}^{-3}$	$6.99\times {10}^{-2}$	1.79
2	4	$1.95\times {10}^{-3}$	$1.76\times {10}^{-2}$	1.99
2	4	$9.77\times {10}^{-4}$	$4.84\times {10}^{-3}$	1.86

Table 3. Table 3: Results Using Third-Order Approximation, without Extra Boundary Refinement ( q b subscript 𝑞 𝑏 q_{b} : number of quadrature points per boundary cell, q v subscript 𝑞 𝑣 q_{v} : number of quadrature points per volume cell, e 0 subscript 𝑒 0 e_{0} : relative error at t = δ ​ t 𝑡 𝛿 𝑡 t=\delta t ).

$q_b$	$q_v$	$\delta x$	$e_0$	EOC
3	9	$1.56\times {10}^{-2}$	$2.11\times {10}^{-1}$	–
3	9	$7.81\times {10}^{-3}$	$3.73\times {10}^{-2}$	2.50
3	9	$3.91\times {10}^{-3}$	$1.04\times {10}^{-2}$	1.84
3	9	$1.95\times {10}^{-3}$	$2.36\times {10}^{-3}$	2.13
3	9	$9.77\times {10}^{-4}$	$5.58\times {10}^{-4}$	2.08

Table 4. Table 4: Results Using Third-Order Approximation, with Extra Boundary Refinement ( q b subscript 𝑞 𝑏 q_{b} : number of quadrature points per boundary cell, q v subscript 𝑞 𝑣 q_{v} : number of quadrature points per volume cell, e 0 subscript 𝑒 0 e_{0} : relative error at t = δ ​ t 𝑡 𝛿 𝑡 t=\delta t ).

$q_b$	$q_v$	$\delta x$	$e_0$	EOC
3	9	$1.56\times {10}^{-2}$	$2.11\times {10}^{-1}$	–
3	9	$7.81\times {10}^{-3}$	$3.45\times {10}^{-2}$	2.61
3	9	$3.91\times {10}^{-3}$	$4.42\times {10}^{-3}$	2.96
3	9	$1.95\times {10}^{-3}$	$6.99\times {10}^{-4}$	2.66
3	9	$9.77\times {10}^{-4}$	$1.04\times {10}^{-4}$	2.75

Table 5. Table 5: Temporal Convergence Tests ( δ ​ t 𝛿 𝑡 \delta t : time step size, e f subscript 𝑒 𝑓 e_{f} : relative error at final time).

$\delta t$	$e_f$	EOC
$0.1$	$2.92\times {10}^{-1}$	–
$0.05$	$1.42\times {10}^{-1}$	$1.04$
$0.025$	$7.53\times {10}^{-2}$	$0.92$

Table 6. Table 6: Condition Numbers of FEM Matrices (2nd Order).

$\delta x$	$N_{DOFs}\left(FEM\right)$	$\kappa \left(FEM\right)$
1.11× 10-1	16	9.79× 103
6.26× 10-2	50	1.24× 104
3.32× 10-2	178	2.15× 104
1.71× 10-2	674	6.82× 104
8.64× 10-3	2626	1.12× 105

Table 7. Table 7: Condition Numbers of IEM Matrices (4th Order).

$\delta x$	$N_{DOFs}\left(IEM\right)$	$\kappa \left(IEM\right)$
1.11× 10-1	56	9.86× 102
6.26× 10-2	100	9.87× 102
3.32× 10-2	200	9.88× 102
1.71× 10-2	364	9.88× 102
8.64× 10-3	720	9.89× 102

Equations

$${\mathcal{E}}_{GL}[\varphi ]={\int }_{\Omega }\frac{ϵ}{2}\Vert \nabla \varphi (\mathbf{x}){\Vert }^2+\frac{[\varphi (\mathbf{x}{)}^2-1{]}^2}{4ϵ}d\mathbf{x},$$

$$\frac{\partial \varphi (\mathbf{x},t)}{\partial t}$$

$$\mu (\mathbf{x},t)$$

$$\varphi (\mathbf{x},0)$$

$$\mu (\mathbf{x},0)$$

$${\mathcal{E}}_S={\int }_{\partial \Omega }\gamma (\varphi (\mathbf{x}))d{s}_{\mathbf{x}},$$

$$\frac{\partial \varphi (\mathbf{x},t)}{\partial t}$$

$${\partial }_n\mu (\mathbf{x},t)$$

$$\frac{{\varphi }^{n+1}(\mathbf{x})-{\varphi }^n(\mathbf{x})}{\delta t}$$

$${\mu }^{n+1}(\mathbf{x})$$

$$\frac{{\varphi }^{n+1}(\mathbf{x})-{\varphi }^n(\mathbf{x})}{\delta t}$$

$${\partial }_n{\mu }^{n+1}(\mathbf{x})$$

$$({\Delta }^2-b\Delta +c){\varphi }^{n+1}(\mathbf{x})$$

$$(\Delta -b){\varphi }^{n+1}(\mathbf{x})+\frac{1}{ϵ}{\mu }^{n+1}(\mathbf{x})$$

$$({\partial }_n+c){\varphi }^{n+1}(\mathbf{x})$$

$$\frac{1}{ϵ}{\partial }_n{\mu }^{n+1}(\mathbf{x})$$

$$f_1(\mathbf{x})$$

$$f_2(\mathbf{x})$$

$$h(\mathbf{x})$$

$$(\Delta -{\lambda }_1^2)(\Delta -{\lambda }_2^2)G_0(\mathbf{x},\mathbf{y})=\delta (\mathbf{x}-\mathbf{y}).$$

$$(\Delta -{\lambda }_i^2)G_i(\mathbf{x},\mathbf{y})=\delta (\mathbf{x}-\mathbf{y})$$

$$G_i(\mathbf{x},\mathbf{y})=-\frac{1}{2\pi }{K}_0({\lambda }_{i}r),i=1,2,$$

$$G_0(\mathbf{x},\mathbf{y})=-\frac{1}{2\pi }\frac{1}{{\lambda }_1^2-{\lambda }_2^2}[K_0({\lambda }_{1}r)-K_0({\lambda }_{2}r)].$$

$$G_0(x,y)=-\frac{1}{4\pi \lambda }r{K}_1(\lambda r),$$

$$G_0(\mathbf{x},\mathbf{y})$$

$$=\frac{1}{{\lambda }_1^2-{\lambda }_2^2}\left[G_1(\mathbf{x},\mathbf{y})-G_2(\mathbf{x},\mathbf{y})\right],$$

$$\Delta G_0(\mathbf{x},\mathbf{y})=\frac{1}{{\lambda }_1^2-{\lambda }_2^2}\left[\Delta G_1(\mathbf{x},\mathbf{y})-\Delta G_2(\mathbf{x},\mathbf{y})\right].$$

$$\Delta G_0(\mathbf{x},\mathbf{y})$$

$$=\frac{{\lambda }_1^2}{{\lambda }_1^2-{\lambda }_2^2}{G}_1(\mathbf{x},\mathbf{y})-\frac{{\lambda }_2^2}{{\lambda }_1^2-{\lambda }_2^2}{G}_2(\mathbf{x},\mathbf{y})$$

$$={\lambda }_1^2{G}_0(\mathbf{x},\mathbf{y})+G_2(\mathbf{x},\mathbf{y})$$

$$={\lambda }_2^2{G}_0(\mathbf{x},\mathbf{y})+G_1(\mathbf{x},\mathbf{y}),$$

$$(\Delta -{\lambda }_1^2)G_0(\mathbf{x},\mathbf{y})=G_2(\mathbf{x},\mathbf{y}),\text{and }(\Delta -{\lambda }_2^2)G_0(\mathbf{x},\mathbf{y})=G_1(\mathbf{x},\mathbf{y}).$$

$$(\Delta -{\lambda }_2^2)(\Delta -{\lambda }_1^2)G_0(\mathbf{x},\mathbf{y})=(\Delta -{\lambda }_1^2)(\Delta -{\lambda }_2^2)G_0(\mathbf{x},\mathbf{y})=\delta (\mathbf{x}-\mathbf{y}).$$

$$S_i[\sigma ](\mathbf{x})={\int }_{\partial \Omega }{G}_i(\mathbf{x},\mathbf{y})\sigma (\mathbf{y})d{s}_{\mathbf{y}}(\mathbf{x}\in {\mathbb{R}}^2,i=0,1,2).$$

$$V_i[f](\mathbf{x})={\int }_{\Omega }{G}_i(\mathbf{x},\mathbf{y})f(\mathbf{y})d\mathbf{y}(\mathbf{x}\in {\mathbb{R}}^2,i=0,1,2).$$

$$\underset{\mathbf{x}\to {\mathbf{p}}^{\pm }}{\lim }\frac{\partial S_i[\sigma ](\mathbf{x})}{\partial {\mathbf{n}}_{\mathbf{p}}}=\left(\pm \frac{1}{2}I+{\partial }_n{S}_i\right)[\sigma ](\mathbf{p})(i=1,2),$$

$${\partial }_n{S}_i[\sigma ](\mathbf{p})=\mathrm{PV}{\int }_{\partial \Omega }\frac{\partial G_i(\mathbf{p},\mathbf{y})}{\partial {\mathbf{n}}_{\mathbf{p}}}\sigma (\mathbf{y})d{s}_{\mathbf{y}}(i=1,2),$$

Full text

1 Introduction

The Cahn-Hilliard equation is frequently used in the phase field model in modeling processes involving the evolution of interfaces ([1, 2, 3]). There has been a large body of work on numerical methods for the Cahn-Hilliard equations (cf. [4, 5, 6, 7] and references therein). Most of the existing methods are finite element/finite difference methods, or spectral methods for simple geometries. Being a fourth order partial differential equation, standard finite difference or finite element methods lead to very ill-conditioned matrices. This in turn causes low accuracy and may require preconditioning for efficient solution. Spectral methods are superior in terms of convergence order and number of unknowns in the discretized system at the cost of only being able to handle simple geometries.

In this paper, we present a high-order accurate numerical method along with an efficient solution algorithm for solving the two-phase Cahn-Hilliard equation using integral equation methods in complex geometry. At each time step, our method requires the evaluation of volume potentials on an adaptive volume grid, followed by the solution of a well-conditioned system of second kind integral equations (SKIEs) with unknown boundary densities. The volume potentials are evaluated in linear complexity with respect to the number of quadrature nodes using the volume fast multipole method (FMM) based on a “box FMM” similar to, for example, [8, 9, 10, 11, 12]. The boundary integral equations involve weakly singular integrals, which are discretized using the quadrature-by-expansion (QBX) method [13]. The idea of QBX is to exploit the smoothness of the (smooth) layer potential away from the surface by forming locally-valid expansions which are then evaluated to compute the near or on-surface value of the potential. Finally, the resulting linear system is solved by iterative solvers such as GMRES [14] with matrix-vector products accelerated by a version of the FMM [15, 16, 17, 18], yielding linear complexity with respect to the number of boundary nodes. When an adaptive volume mesh is used, the adaptive box structure is adjusted after each step to capture the moving interface by refining in the neighborhood of the interface and coarsening away from the interface based on solution gradient; in addition, local refinement near the domain boundary is used to achieve higher order accuracy during volume potential evaluation. The overall scheme is well conditioned, potentially high-order in space, and has asymptotically optimal complexity in terms of problem size.

The outline of the remainder of the paper is as follows. In Section 2, we present the model equations and the temporal discretization scheme. In Section 3, we analyze the static problem at each time step, and provide a mathematical basis for the SKIE formulation. In Section 4, we present the SKIE formulation for the static problem. In Section 5, we present the numerical algorithms to accompany the method formulation so that the final solver is robust and has linear complexity, and in Section 6 we present some numerical results. In Section 7, we summarize the paper and discuss some future perspectives.

2 The Cahn-Hilliard Equation: Original Problem Specification

We consider a two-phase Cahn-Hilliard equation derived from a simple phase-field model given by the following Ginzburg-Landau energy functional [19]:

[TABLE]

where $\phi(\mathbf{x})$ is the composition field and $\epsilon$ represents the interface thickness. When $\Omega=\mathbb{R}^{2}$, taking the gradient flow of (1) in $H^{-1}(\mathbb{R}^{2})$ yields the Cahn-Hilliard equation [20]:

[TABLE]

where $\mu$ is the Fréchet derivative of (1), also known as the chemical potential.

In practice, we are interested in solving the initial-boundary value problem on bounded domains. The initial conditions are

[TABLE]

Being a fourth order PDE, the Cahn-Hilliard equation requires two sets of boundary conditions. In the wetting problem, the boundary models solid surfaces, and a physically relevant way to impose boundary conditions is by adding a surface energy term to the free energy to account for the interaction with the solid wall: [21]

[TABLE]

where we choose $\gamma(\phi)=\frac{\sqrt{2}}{3}\cos\theta_{Y}\sin\left(\frac{\pi}{2}\phi\right)$, which gives the equilibrium contact angle (Young’s angle) of $\theta_{Y}$ [22]. Re-take the gradient flow but now with the total free energy $\mathcal{F}=\mathcal{E}_{GL}+\mathcal{E}_{S}$, and the extra boundary energy term yields the relaxation boundary condition [23]

[TABLE]

as well as the zero-flux boundary condition

[TABLE]

where $\partial_{n}=\mathbf{n}\cdot\nabla$ with $\mathbf{n}$ being the unit outward normal vector of $\partial\Omega$.

Using Rothe’s method, an energy decaying time stepping scheme for the system (2 – 5) and (7) – (8) is introduced in [24], which is based on the idea of convex splitting [25]: at the n-th time step, given $\phi^{n}(\mathbf{x})$ , find $\phi^{n+1}(\mathbf{x})$ such that for all $\mathbf{x}$ inside the domain,

[TABLE]

and that for all $\mathbf{x}$ on the boundary,

[TABLE]

The splitting parameter $s$ is chosen to be large enough such that the scheme is energy-stable (under the assumption that $\|\phi\|_{\infty}<\infty$, see [24] for details). If a bound to the solution is known $\|\phi\|_{\infty,\Omega}\leq M$, $\forall t>0$, then letting $s\geq(3M^{2}-1)/2$ guarantees energy stability of the scheme; however, even though the numerical solutions mostly stay around the $[-1,1]$ range, the mathematical problem of solution boundedness for Cahn-Hilliard equation is still open. In this paper we choose to fix $s=1.5$, which is sufficient to allow for arbitrary $\delta t$ in all our tests.

It is worth noting that time discretization of the Cahn-Hilliard equation is still an active research area. In this paper, we choose to adopt the first order linearly-implicit convex-splitting method and direct our main efforts to solving the static forced modified biharmonic subproblems using integral equation approach. Other time discretization methods include: higher-order convex splitting methods ([26, 27]) , stabilization methods ([7, 28]) , the method of invariant energy quadratization (IEQ) ([29, 30]), and the recent scalar auxiliary variable (SAV) method [31]. Some of those advanced time integration schemes (like the SAV method) produce fourth-order static boundary value problems with constant coefficients, similar to the one treated in this paper. The integral equation formulation and fast algorithms in this paper can thus be straightforwardly generalized to adopt such time discretization techniques.

For convenience of deriving integral representations of the solution, we rewrite the problem (9)-(12) by collecting unknowns onto one side of the equation:

[TABLE]

where $b=\frac{s}{\epsilon^{2}},c=\frac{1}{\epsilon\delta t}$, and the inhomogeneous terms $f_{1},f_{2}$ and the boundary data $h$ are given by

[TABLE]

3 Fundamental Solutions and Their Properties

Denote the roots of the quadratic equation $x^{2}-bx+c=0$ by $\lambda_{1}^{2}$ and $\lambda_{2}^{2}$. Then the fourth order operator $\Delta^{2}-b\Delta+c$ in (13) can be factored as $(\Delta-\lambda_{1}^{2})(\Delta-\lambda_{2}^{2})$. To employ an IE method for the solution of (13), we seek a fundamental solution $G_{0}$ that should solve, in a weak sense,

[TABLE]

Motivated by the factorization structure, we recall the Green’s function of the Yukawa operator that satisfies (weakly)

[TABLE]

and is given by the expression

[TABLE]

where $r=\|x-y\|_{2}$ and $K_{0}$ is the modified Bessel function of the second kind of order zero (see, for example, [32]). For (20), we formulate the fundamental solution as follows:

Definition 1**.**

The fundamental solution $G_{0}$ is given by the formula

[TABLE]

In the special case of $\lambda_{1}=\lambda_{2}=\lambda$, the funcdamental solution $G_{0}$ is defined by taking the limit $\lambda_{1}\rightarrow\lambda_{2}$ of (23), and is given by the formula

[TABLE]

where $K_{1}$ is the modified Bessel function of the second kind of order one.

For simplicity of discussion, we assume that $\lambda_{1}^{2}\neq\lambda_{2}^{2}$ for the rest for this paper. In practice, if $\lambda_{1}^{2}=\lambda_{2}^{2}$, the equality can always be broken by slightly changing the values of the time discretization parameters $s$ and $\delta t$.

Proposition 2**.**

$G_{0}$* as defined by (23) satisfies (20).*

Proof.

By definition (23),

[TABLE]

therefore,

[TABLE]

Using (21), we have

[TABLE]

that is,

[TABLE]

Combining (25) and (21), we have

[TABLE]

∎

Now we are equipped to define layer and volume potential operators.

Definition 3** (Single layer potentials).**

Given a density function $\sigma(\mathbf{x})\in C(\partial\Omega)$, the single layer potential operators are

[TABLE]

Definition 4** (Volume potentials).**

Given a density function $f(\mathbf{x})\in C(\Omega)$ , the volume potential operators are

[TABLE]

Note that here we only require the density functions to be continuous. In practice, higher regularity is desired for use of high order discretization, and our solver in this paper builds on $C^{1}$ density functions.

Regarding the jump conditions of layer potentials across the domain boundary, we have the following result:

Theorem 5** (Jump relations).**

The normal derivative of $S_{0}[\sigma]$ has no jump across the boundary $\partial\Omega$. However, when $\mathbf{x}$ approaches a point $\mathbf{p}\in\partial\Omega$ nontangentially, the normal derivative of the single layer potential operators $S_{i}[\sigma]$ ($i=1,2$) on $\partial\Omega$ satisfies the following jump relation:

[TABLE]

where $\mathbf{n}_{\mathbf{p}}$ is the unit outward normal of $\partial\Omega$ at $\mathbf{p}$, and

[TABLE]

and $I$ is the identity operator. Here $\mathbf{x}\rightarrow\mathbf{p}^{\pm}$ means that $\mathbf{x}$ approaches $\mathbf{p}$ from the exterior $(+)$ or the interior $(-)$ of the domain, respectively, and $\mathop{\mathrm{P}V}\cdot$ denotes the Cauchy principal value.

Proof.

We expand $K_{0}$ at the origin using asymptotic formulae [32, (9.6.10)–(9.6.13)]. When $z\rightarrow 0$,

[TABLE]

where $\psi(x)$ is the psi (digamma) function defined as

[TABLE]

and $\gamma=0.5772156649\ldots$ is Euler’s constant [32].

From (30) it is clear that the leading-order non-smooth behavior results from the logarithmic term. Then (28) follows directly from the well-known jump condition for the normal derivative of the single layer potential with logarithmic kernel ([33, Lemma 3.30]). ∎

Using the compactness of integral operators with weakly singular kernels [34, Theorem 1.11], we also have the following corollary regarding the compactness of those single layer potential operators and their derivatives

Corollary 6** (Compactness).**

For a Lipschitz domain $\Omega$, all three single layer potential operators are compact in $C(\Omega)$. Also, the principal value parts of $\partial_{n}S_{i}[\sigma]$ ($i=0,1,2$) are all compact.

4 Integral Equation Formulation

In this section, we present a second-kind integral equation formulation for the system (13)–(19). First we introduce a construction of the volume potential that reduces the problem to a pure (volume-homogeneous) boundary value problem (PBVP). Then we present an integral equation formulation for the PBVP that only requires solving a Fredholm integral equation of the second kind.

4.1 Reduction to Pure Boundary Value Problem

We choose to represent the solutions $\phi^{n}$, $\mu^{n}$ in terms of sums of volume and layer potentials. We first introduce the volume potential in the terms of $\tilde{u}$ and $\tilde{v}$:

[TABLE]

where we require that $V[\phi^{n}]$ satisfies

[TABLE]

and $V[\cdot]$ is a volume potential operator to be defined.

For the sole purpose of removing inhomogeneities in (13) and (14), one could simply use the volume potential

[TABLE]

In numerical computation, however, evaluating $f_{1}$ (cf. (17)) directly is not desirable since it involves second order differentiation and thus requires $\phi^{n}\in C^{2}$, which can be overly restrictive. Applying Green’s second identity to (35) yields

[TABLE]

Our construction of $V[\phi^{n}]$ from (32) is obtained by dropping the last two boundary integral terms and making use of (25), leading us to define

[TABLE]

Inserting (36) into (34) shows that the volume condition is satisfied despite the dropping of the boundary terms. By using our construction (36) numerical differentiation is completely avoided in the volume potential evaluation stage.

Clearly, one does not need to solve for but only evaluate $\tilde{u}$ and $\tilde{v}$. Specifically, given $\phi^{n}$, we first evaluate $\tilde{u}$ by computing the volume integral $V[\phi^{n}]$ directly with the help of the box FMM. In the same FMM pass, we also obtain $\Delta\tilde{u}$ by taking derivatives of the local expansions. The loss of accuracy incurred due to numerical differentiation can be recovered by increasing the FMM order. Then we can compute $\tilde{v}$ directly from its defining formula (33).

Then $u$ and $v$ are the solutions of the following pure boundary value problem

[TABLE]

where the boundary data $g_{1}$ and $g_{2}$ are given by

[TABLE]

4.2 Second Kind Integral Equation Formulation for the Pure Boundary Value Problem

We now derive a second kind integral equation (SKIE) formulation for the boundary value problem 37 to 42. We first represent $u$ by the formula

[TABLE]

where $\sigma_{i}$ ($i=1,2$) are unknown densities on $\partial\Omega$. Obviously, this representation satisfies (37). Substituting (43) into (38) yields

[TABLE]

Note that since $b=\lambda_{1}^{2}+\lambda_{2}^{2}$, the following corollary follows naturally from the definition of the Yukawa kernels (21):

Corollary 7**.**

When $\mathbf{y}\neq\mathbf{x}$, $G_{i}(\mathbf{x},\mathbf{y})$ ($i=1,2$) are both locally eigenfunctions of the Yukawa operator $\Delta-b$. Specifically,

[TABLE]

Proof.

We only prove the first identity. The other one follows analogously. Substituting $b=\lambda_{1}^{2}+\lambda_{2}^{2}$ into the left hand side of the identity,

[TABLE]

By definition, $(\Delta-\lambda_{1}^{2})G_{1}(\mathbf{x},\mathbf{y})=\delta(\mathbf{x},\mathbf{y})=0$. Thus

[TABLE]

∎

On the other hand, (25) yields

[TABLE]

Thus

[TABLE]

Combining the jump relations (28), the boundary conditions (39) and (40), and the representations (43) and (47), we obtain the following system of boundary integral equations:

[TABLE]

where

[TABLE]

and the entries of the operator matrix $A$ are given by the formulae

[TABLE]

Clearly, $D$ has nonzero determinant, and all entries of $A$ are compact operators (Corollary 6). Therefore, the system (48) is a Fredholm integral equation of the second kind. When discretizing a second kind integral equation, the condition number of the resulting linear system remains bounded when the number of unknowns increases as the discretization gets refined. For an iterative solver achieving a fixed residual norm, the numerical solution thus remains accurate as the mesh is refined. Besides, with bounded condition number, the number of GMRES iterations is less likely to grow when refining the mesh.

5 Numerical Realization: Discretization and Algorithms

For the implementation of our method, the volume potentials in (36) are evaluated using a version of the the volume FMM ([8, 9, 10, 11, 12]). We further discretize the SKIE system (48) using the Nyström method using QBX quadrature and solve the linear system with the help of GMRES. The required matrix-vector products are carried out using the GIGAQBX quadrature-enabled fast algorithm [13, 35, 36]. To handle complex input geometries, we use Gmsh [37] to generate a triangulation along the domain boundary with each panel roughly of the size $h_{b}$. Over each boundary panel, Legendre-Gauss quadrature nodes of order $q_{b}$ are used to form the discretization.

If $\Omega=[x_{l},x_{r}]^{2}$, the volume FMM bridges seamlessly with the boundary discretization; however, additional treatment is needed when dealing with general geometries, as the volume FMM requires box-shaped partitions to operate. In this section, we present those extra steps necessary to achieve flexibility in handling complex geometries.

5.1 Regularity and Non-Square Domains

For non-square domains, the volume potentials are evaluated over a box $B$ that encloses the physical domain such that $\Omega\subset B$. Use of the box FMM requires that volume source densities be extended to $B\setminus\Omega$. Noting that $f_{1},f_{2}$ are all algebraic functions of $\phi^{n}$, we only need to extend $\phi$ to $B\setminus\Omega$. There are several factors that need to be considered when deciding how to extend $\phi$:

The smoother the extension is, the higher accuracy in the order of approximation for those densities can be achieved using (piecewise) polynomials, yielding better overall accuracy.

• 2.

Higher order extensions are more costly to evaluate. Typically for $C^{k}$ extension one needs to solve a PDE where $k+1$ boundary conditions can be imposed, e.g., $\Delta^{k+1}u=0$. It is natural to ask that the extension process does not cost much more than solving the original problem, i.e., $k\leq 1$.

• 3.

Our SKIE formulation poses certain minimum regularity requirements. Specifically, evaluation of $g_{2}$ in (40) requires at minimum $\tilde{v}\in C^{1}(B)$, which then requires $\tilde{u}$ to be in $C^{3}(B)$.

To satisfy the condition $\tilde{u}\in C^{3}(B)$, the density $\phi$ should satisfy certain regularity requirements. First we note the standard regularity result for the Poisson equation [38]:

Theorem 8**.**

Let $\Omega\subset\mathbb{R}^{d}$ be open and bounded,

[TABLE]

where $\Phi$ is the fundamental solution. Then if $f\in C_{0}^{\alpha}(\overline{\Omega})$, $0<\alpha<1$, then $u\in C^{2,\alpha}(\overline{\Omega})$, and

[TABLE]

The leading order singular (least regular) term in our $V[\phi^{n}]$ construction is $V_{1}[f_{2}]$, for which we have the following result:

Proposition 9**.**

If $f_{2}\in C^{1}(B)$, the volume potential $V_{1}[f_{2}]\in C^{3}(B)$.

Proof.

Since, like the Poisson kernel, the kernel $G_{1}$ also has a leading-order logarithmic singularity, Theorem 8 also applies to $V_{1}[f_{2}]$. Combined with the assumption that the source density is in $C^{1}(B)$, the proposition follows. ∎

Based on these considerations, we choose to use $C^{1}$ extensions. To obtain such an extension, we assume that $\Omega$ is simply connected and solve the following Stokes problem:

[TABLE]

where $\mathbf{u}$ is the velocity field of the fluid and $p$ is the pressure field. The Stokes problem (53, 54, 55) can also be solved via an SKIE formulation [39]. We evaluate the stream function $\omega$ up to a constant by evaluating complex layer potential representations in [40], such that

[TABLE]

then $\omega$ is a biharmonic function that solves the boundary value problem

[TABLE]

Since $\partial\Omega$ is simply connected, we only need to add a constant to $\omega$ to obtain a $C^{1}$ extension of $\phi$:

[TABLE]

Remark 10**.**

The assumption of $\Omega$ being simply connected is for simplicity of the formulation and not critical. If $\Omega$ is topologically more complicated, methods in [41, 40] can still be used for performing $C^{1}$ extension.

5.2 Boundary Layers

From matched asymptotic analysis of the Cahn-Hilliard equation (2)–(8) [42], it is apparent that for very small $\epsilon$, the solution can develop boundary layers.

Boundary layers cause stability concerns for our scheme using $C^{1}$ extension. Since the exterior biharmonic problem does not admit a Laplace-type maximum principle, the presence of the boundary layer causes the extended density to be very large in magnitude compared to the density inside the domain, resulting in lost accuracy after each time step. As a result, the numerical solution may not remain bounded after several time steps and lose energy stability.

Dealing with such boundary layers numerically is the subject of ongoing research. Using an adaptive volume mesh that refines towards the boundary can be helpful. Instead of evaluating the volume potentials over one single bounding box $B$, we can choose $B$ to be union of a set of boxes that covers $\Omega$. By reducing the Hausdorff distance between $\partial B$ and $\partial\Omega$, the extended density is less likely to blow up. However, such refinement is very expensive in practice, and there is no guarantee that the boundary layer is resolved a-priori since a reliable refinement criterion is not known.

For long-time simulations, we provide an alternative formulation which is a less costly approximation based on the total energy of the system. When solving the SKIE (48), we use a homogeneous right hand side for the second equation $g_{2}=0$, which effectively makes the layer potential $u$ satisfy the contact angle dynamics, while ignoring the second boundary condition. After solving for $\sigma_{1}$ and evaluating $u$, we add a constant to $u$ such that the integral of $\phi$ over the whole domain is conserved. This approach conserves the total mass and keeps the numerical solution bounded at all time.

Formally, the approximation replaces the representation (43) with the following stabilized representation:

[TABLE]

where $C$ is the added constant. The net effect of the approximation spreads the influence of boundary layers across the whole domain. Long-time simulations using the stabilized representation are performed and the results are in line with those obtained from finite element method. Although the modification seems crude, it maintains the most relevant solution characteristics including:

$H^{-1}$ gradient flow structure in the bulk region, since the volume potential is unchanged.

• 2.

Contact angle dynamics from the first boundary condition.

• 3.

Integral conservation over the whole domain.

6 Numerical Results

6.1 Spatial Convergence

Based on error analysis of the volume FMM ([9]), QBX ([43]), the GIGAQBX fast algorithm ([13], [36]) and standard approximation theory ([44]), we expect the following error estimate for the overall numerical accuracy of our scheme:

Heuristic 1** (Numerical Accuracy).**

Assuming density extension has regularity $\phi^{0}\in\mathcal{C}^{k}(B)$, denote the numerical solution at $t=\delta t$ by $\phi_{h}^{1}$ and its relative error by $E_{0}:=\lVert\phi_{h}^{1}-\phi(\delta t)\rVert_{\infty}/\lVert\phi(\delta t)\rVert_{\infty}$. Asymptotically, for a single time step, we have

[TABLE]

where $C_{0}$ is a constant independent of $\delta_{t}$, $E_{v}$ includes the errors from volume potential evaluation and $E_{b}$ includes the errors from the BIE solution.

For $E_{v}$, we have

[TABLE]

where $h_{v}$ is the volume mesh size, $q_{v}$ is the number of quadrature points per volume cell, $m_{v}$ is the volume FMM order, $\xi$ is the error of near-field volume potentials (controlled by the tolerance of adaptive quadrature over nearfield boxes), and $C_{1},C_{2}$ are constants independent of $h_{v}$ and $\xi$. 2. 2.

For $E_{b}$, assuming that the solution density $\sigma$ is smooth and that its Hölder norms $\lVert\sigma\rVert_{C^{s,\beta}}$ for $s\in\mathbb{N},\beta\in[0,1)$ can be controlled with $\lVert\phi(\delta t)\rVert_{\infty}$, we have

[TABLE]

where $h_{b}$ is the boundary mesh size, $q_{b}$ is the number of quadrature points per boundary cell, $p$ is the QBX order, $m_{b}$ is the FMM order used in GIGAQBX, $\eta$ is the error of the linear solver (controlled by the tolerance of GMRES), and $C_{3},C_{4},C_{5},C_{6},C_{7}$ are constants independent of $h_{b}$ and $\eta$.

Among all the error sources listed in Heuristic 1, we are most interested in the discretization errors that decay algebraically with mesh refinement. In this section, we describe experiments to verify this expected error decay, by choosing other parameters of the scheme so as to minimize the contributions of their corresponding errors. Specifically, we fix the tolerance for GMRES to be $10^{-14}$, and the tolerance for (precomputed) near-field adaptive quadrature in the volume FMM to be $10^{-13}$. Also, we set the FMM orders and the QBX order to be sufficiently high. Specifically, we let the QBX order match the boundary quadrature order $q_{b}=p+1$, and fix the volume FMM order and the GIGAQBX’s FMM order $m_{v}=m_{b}=8$. We use two-dimensional Taylor expansions for the FMM, thus achieving FMM order $8$ in $2$D requires $\frac{1}{2}\times 8\times 9=36$ terms. To back our choices, neither increasing the QBX order $p$ up to $10$ nor the FMM order $m_{v},m_{b}$ up to $10$ yields obvious effects on the results, indicating that the error in the computation of the volume potential is dominated by discretization errors. The same characteristic mesh size $h_{b}=h_{v}=\delta x$ is used for volume and boundary meshes. For simplicity, we also keep the orders of volume and boundary discretizations the same.

To gauge convergence, we fix $\delta t=10^{-2}$ and perform self-convergence tests by refining both spatial grids successively and calculating the relative error from the result of one time step $e_{0}:=\lVert\phi^{1}_{h}-\phi^{1}_{0}\rVert_{\infty}/\lVert\phi^{1}_{h}\rVert_{\infty}$, where $\phi^{1}_{0}$ is the reference solution. The tests are performed on a circular domain with radius $0.247$ centered at the origin. The bounding box is chosen to be $B=[-0.25,0.25]^{2}$. The initial condition is

[TABLE]

where $\epsilon=10^{-2}$, and $\mathbf{e}_{1}=(1,0)$.

In this work, since we are using $C^{1}$ extension of $\phi$, based on the estimates above, when the discretization error dominates, the total error of the scheme should scale as $\mathcal{O}(h_{b}^{q_{b}})+\mathcal{O}(h_{v}^{\min(\sqrt{q_{v}},2)})$ when $h_{b},h_{v}\rightarrow 0$. In particular, when $h_{b}=h_{v}=\delta x$, and $q_{b}=\sqrt{q_{v}}=r$, the total discretization error is expected to be $\mathcal{O}(\delta x^{r})$ as $\delta x\rightarrow 0$.

When using piecewise constant approximations for both boundary and volume discretizations, the results are shown in Table 1, where $q_{b},q_{v}$ stand for the number of quadrature points per cell on the boundary mesh and the volume mesh, respectively. The convergence order is calculated using values of the error measure $e_{0}$ from the current row and the row above. The same statistics for the second order discretizations (i.e. piecewise linear approximations on the boundary mesh and piecewise bi-linear approximations on the volume mesh) are shown in Table 2.

In both cases, the numerical results align well with expectations. When $r=1$, $q_{b}=q_{v}=1$, discretizations for both boundary and volume densities are piecewise constant (first order), and Table 1 demonstrates empirical convergence order close to $1$. When $r=2$, $q_{b}=2,q_{v}=4$, boundary densities are discretized with piecewise linear functions, and volume densities are discretized with piecewise bilinear functions, in which case Table 2 shows empirical convergence order close to $2$.

When higher-order volume approximations are used with a uniform volume mesh, the convergence order will be limited to $2$ due to the volume densities being limited to only $C^{1}$ at the boundary. For example, as shown in Table 3.

However, higher order approximations can still be used by adding extra levels of refinement to volume cells that intersect $\partial\Omega$. To guide the mesh refinement process, we form a two-dimensional tensor-product Legendre series of degree $q_{v}$ (using the discrete tensor-product Legendre transform) over each cell, and compute the sum of the absolute values of the $q_{v}+1$ coefficients of the leading terms (terms of total degree $q_{v}$) as the error indicator for that cell. To make the adaptive mesh for the box FMM, we start from a uniform mesh with prescribed $\delta x$, and iteratively perform refinement cycles as follows:

Compute the average values of the error indicators for the interior cells as $M_{i}$, and that for the boundary-intersecting cells as $M_{b}$. 2. 2.

Refine each boundary-intersecting cell by one level (into four cells). 3. 3.

Refine some additional cells as necessary to keep the tree balanced.

We repeat the refinement cycles until $M_{b}\leq M_{i}$. As shown by the results in Table 4, since the source density is smooth away from the boundary, by performing such refinements, the solver is able to take advantage of higher order of convergence in the bulk region where the source density is smooth.

6.2 Temporal Convergence

In terms of temporal convergence, we expect to have first order in time due to use of the convex splitting scheme of [24]. We still use the initial condition (61) over a disk domain to do perform this test. To verify the expectation, we fix $\epsilon=0.5$ and the final time $t_{f}=0.1$, and perform self-convergence tests by taking different $\delta t$. The numerical error is measured at $t=t_{f}$ using the formula $e_{f}:=\lVert\phi^{f}_{h}-\phi^{f}_{0}\rVert_{\infty}/\lVert\phi^{f}_{h}\rVert_{\infty}$, against the reference solution obtained by taking $\delta t=0.0125$. The results are presented in Table 5, confirming that the scheme is of first order in time.

6.3 Short-time Dynamics

In this example, we test for the ability of our scheme to resolve small structures using high order approximation in combination with adaptive mesh refinement. We mimic, using the initial condition, a near-uniform state with small random perturbations. Those perturbations will lead to a phase separation process. We expect that our scheme is able to resolve the dynamics of structures at all scales during the phase separation.

The setup tries to capture short-time dynamics of the Cahn-Hilliard equation. We choose $\epsilon=10^{-2}$, $\delta t=10^{-4}$, and use fourth-order quadrature rules for volume and boundary discretizations. For volume densities,

[TABLE]

(where $P_{k}$ stands for Legendre polynomial of degree $k$ defined on the concerning cell) we use

[TABLE]

as the error indicator for the cell.

The simulation starts with an initial profile generated using the following procedure:

A uniform $200\times 200$ Cartesian grid is generated on $B$. Then for each grid point a real number is drawn from a uniform random distribution on $[-10^{-3},10^{-3}]$. 2. 2.

A cubic B-spline representation of the initial $\phi$ is formed using the Cartesian grid points as knots and the random numbers as interpolating values. 3. 3.

The values and derivatives of the intial $\phi$ are then evaluated using the cubic B-spline representation on the quadrature points of the box mesh as well as the boundary mesh.

Numerical solutions are shown in Figure 1. The first $3$ time steps are performed solving the full SKIE for the boundary integral equation to capture the fast dynamics near the boundary, while the later steps are performed using the approximation (60). In our experience, the numerical solution tends to form boundary layers in regions away from contact lines (where its tangential derivative is close to [math]). The color maps are shown for each sub-figure individually, while the contours are all showing the zero level set of $\phi$ over the whole bounding box. The short-time dynamics can be understood as interplay of three mechanisms: surface-energy driven phase transition, bulk-energy driven spinodal decomposition, and solid boundary induced nucleation. In later time steps, the smoothing term in the free energy starts to play a more important role and leads to coarsening of the patterns. As a benefit of the high-order accuracy of our scheme, our method is able to resolve the interfaces as well as the small patterns formed on the boundary very well.

6.4 Long-time Dynamics

Since our solver is only first-order in time, choosing larger time step sizes reduces solution accuracy significantly; however, through this example we demonstrate that even with large time steps, higher spatial order is still useful for preserving certain qualitative properties of the solution. In this example, we verify that the solver can preserve solution symmetry across long time scales.

We still use $\epsilon=10^{-2}$, but take larger time step size $\delta t=0.5$. The initial condition is

[TABLE]

where $L=0.5$ is the size of the bounding box. Note that $\phi^{0}(x,-y)=\phi^{0}(x,y)$. For this example we use (60) from the start. Since $\delta t$ is rather large compared to $\delta x$, the temporal truncation error of convex splitting scheme dominates; however, we expect that the scheme is still stable and reaches an equilibrium state when $t\rightarrow\infty$, and that the solution $\phi$ at any time stays symmetric with respect to the $x$ axis.

Numerical results are shown in Figure 2. The simulation stops at $t=780$ when $\lVert\phi^{n+1}-\phi^{n}\rVert_{\infty}<10^{-14}$, which confirms that the numerical solution reaches an equilibrium state. Besides being able to resolve the interface and the contact lines, our method is also shown to be able to properly resolve the interactions between the contact lines and the geometry with no alignment requirements between the two sets of meshes: in the case circular domain, as long as both the boundary mesh and the box mesh are symmetric, the solution symmetry will be well-preserved from the beginning to the equilibrium state.

Then, we demonstrate the ability of our scheme to handle moderately complex smooth geometries. We report the numerical solution at $t=1$ for the same problem setup, but with different physical domains. The results are shown in Figure 3.

The total energy of the numerical solution for different geometries is shown in Figure 4. According to the theory, the energy decaying property is preserved as long as the numerical solution remains bounded. With the stabilization applied, we observe that the numerical solution remains bounded for various geometries.

6.5 Comparisons with FEM

It is a natural question to ask how the integral equation approach that we take in this article compares to more popular finite element methods. In this example we aim to provide some perspective on this by presenting a comparison between our solver and a simple finite element solver.

The finite element solver we consider is based on a mixed formulation for the time-discretized system (9)–(12), using $Q^{1}$ elements for both $\phi$ and $\mu$. We report the condition number and number of degrees of freedom with their dependence on $\delta x$ for FEM in Table 6, where $\delta x$ is calculated as the square root of average cell measure, and the same statistics for IEM in Table 7. The tests are performed over a circular domain with radius $0.247$. In all tests with both the FEM and the IEM, we fix $\epsilon=0.5$, $\theta_{Y}=\pi/3$ and $\delta t=1$. The rest of the parameters for the IEM are the same as in 6.1. The data confirms that, being a second-kind integral equation, the linear system has bounded condition number that does not scale with mesh size, unlike the FEM whose linear system’s condition number grows in an unbounded manner when refining the mesh. In addition, when $\delta x\rightarrow 0$, the number of unknowns in the linear solve scales as $\mathcal{O}(1/\delta x)$ for our method, while it scales as $\mathcal{O}(1/\delta x^{2})$ for the FEM.

7 Conclusions

We have introduced an integral equation method for the Cahn-Hilliard equation in bounded 2D domains with solid boundary conditions. For each time step, the method consists of two stages: (1) evaluate volume potentials and (2) solve remaining SKIE. Fast multipole based fast algorithms are used to attain linear complexity for both stages. To handle complex geometry, a $C^{1}$ extension of $\phi$ onto a larger domain is performed by solving an additional SKIE. Through numerical experiments, the method is shown to require fewer degrees of freedom and maintain bounded condition number compared to finite element method. To remedy some of the numerical stability concerns for long-time simulations caused by boundary layers, we have proposed a stabilized representation approximating the original problem, and we have validated the approximation with numerical results.

The present work opens some exciting future perspectives. We will seek to apply the same techniques in this paper to second-order time discretization schemes, for example the SAV method [31] offering higher accuracy order in time as well as superior energy stability. In addition, the challenges posed by boundary layers are by no means unique to Cahn-Hilliard equations, and it is of great interest to develop fast volume potential evaluators that remain stable and accurate when the source density develops boundary layers. Once the boundary layer issue can be handled, our SKIE formulation will be immediately applicable to the full long-time simulation.

8 Acknowledgments

This research was supported in part by the National Science Foundation under grants DMS-1418961 and DMS-1654756, and by the Hong Kong Research Grant Council (RGC-GRF grants 605513 and 16302715, RGC-CRF grant C6004-14G, and NSFC-RGC joint research grant N-HKUST620/15). Any opinions, findings, and conclusions, or recommendations expressed in this article are those of the authors and do not necessarily reflect the views of the National Science Foundation or the Hong Kong Research Grant Council; Neither NSF nor HKRGC has approved or endorsed its content. Part of the work was performed while the authors were participating in the HKUST-ICERM workshop ‘Integral Equation Methods, Fast Algorithms and Their Applications to Fluid Dynamics and Materials Science’ held in 2017.