# A CONSTRAINED CONVEX SPLITTING SCHEME FOR THE VECTOR-VALUED CAHN-HILLIARD EQUATION

• LEE, HYUN GEUN (DEPARTMENT OF MATHEMATICS, KWANGWOON UNIVERSITY) ;
• LEE, JUNE-YUB (DEPARTMENT OF MATHEMATICS, EWHA WOMANS UNIVERSITY) ;
• SHIN, JAEMIN (INSTITUTE OF MATHEMATICAL SCIENCES, EWHA WOMANS UNIVERSITY)
• Accepted : 2019.03.19
• Published : 2019.03.25

#### Abstract

In contrast to the well-developed convex splitting schemes for gradient flows of two-component system, there were few efforts on applying the convex splitting idea to gradient flows of multi-component system, such as the vector-valued Cahn-Hilliard (vCH) equation. In the case of the vCH equation, one need to consider not only the convex splitting idea but also a specific method to manage the partition of unity constraint to design an unconditionally energy stable scheme. In this paper, we propose a constrained Convex Splitting (cCS) scheme for the vCH equation, which is based on a convex splitting of the energy functional for the vCH equation under the constraint. We show analytically that the cCS scheme is mass conserving and unconditionally uniquely solvable. And it satisfies the constraint at the next time level for any time step thus is unconditionally energy stable. Numerical experiments are presented demonstrating the accuracy, energy stability, and efficiency of the proposed cCS scheme.

# 1. INTRODUCTION

The Cahn–Hilliard (CH) equation was originally introduced as a phenomenological model of phase separation in a binary alloy [8] and has been applied to a wide range of problems [9]. The CH equation is derived from the Ginzburg–Landau (GL) energy functional:

E(c) := $\int_{\Omega}\left(F(c)+\frac{\epsilon^{2}}{2}|\nabla c|^{2}\right)$ dx,       (1.1)

where Ω is a domain in ℝd (d = 1, 2, 3), c is the concentration field, F(c) is the free energy density for c, and ϵ > 0 is the gradient energy coefficient. The CH equation is a gradient flow for the GL energy functional (1.1) in the H−1 inner product thus the GL energy functional is nonincreasing in time. Since the CH equation cannot be solved analytically in general, numerical methods are commonly used to study the dynamics of the CH equation. Among them, the convex splitting idea [2, 13, 14] has attracted considerable attention, in which the GL energy functional is split into convex and concave parts:

E(c) = Ec(c) − Ee(c) = $\int_{\Omega}\left(F_{c}(c)+\frac{\epsilon^{2}}{2}|\nabla c|^{2}\right) d \mathbf{x}-\int_{\Omega}$ Fe(c) dx,

where Fc(c) and −Fe(c) are convex and concave parts of F(c), respectively. And Ec(c) and Ee(c) are treated implicitly and explicitly, respectively. This idea has been applied to develop uniquely solvable and unconditionally energy stable schemes for a wide class of gradient flows [17, 24, 28, 29, 30, 36, 37, 38].

In order to model phase separation in multi-component systems, several generalizations of the GL energy functional have been introduced and studied [5, 6, 7, 10, 11, 12, 15, 16, 18, 19, 20, 21, 22, 23, 25, 26, 27, 31, 32, 33, 34, 35]. However, there were few efforts on applying the convex splitting idea to multi-component systems. There are two noteworthy related works [21, 32]. In [21], a semi-implicit scheme partially using the convex splitting idea was presented to solve multi-component systems. The semi-implicit scheme in [32] treats all the nonlinear term of multi-component systems explicitly and adds an extra stabilizing term. It can be considered as applying the convex splitting idea when the magnitude of the extra stabilizing term is sufficiently large thus the energy stability depends on the magnitude of the extra stabilizing term.

Unlike the CH equation, designing an unconditionally energy stable scheme for multicomponent systems requires not only the convex splitting idea but also a specific method to manage the partition of unity constraint (the sum of concentration fields must be unity). Note that the importance of the constraint becomes obvious in proving theorem 3.6 in section 3 and some of numerical schemes add the constraint as a part of explicit equations to be solved not as a consequence of the numerical scheme. The authors in [21] solved multi-component systems only for c1, , cN−1, where ci is the concentration field of the phase i and N is the number of phases, and enforced cN = 1−$\sum_{i=1}^{N-1}$cto satisfy the constraint at the next time level. The author in [32] used the Schur complement method to solve multi-component systems for c1, . . . , cN−1 and the constraint c1 + · · · + cN = 1 together.

In this paper, we propose a constrained Convex Splitting (cCS) scheme for the multi– component system used in [5, 10, 12, 18, 19, 20, 21, 22, 23, 27, 31, 32, 33, 34], referred to as the vector-valued Cahn–Hilliard (vCH) equation. This scheme is based on a convex splitting of the energy functional for the vCH equation under the constraint. We show analytically that the cCS scheme is mass conserving and unconditionally uniquely solvable. And it satisfies the constraint at the next time level for any time step thus is unconditionally energy stable. The proposed cCS scheme is a first attempt to achieve unconditional energy stability by applying the convex splitting idea to multi-component systems.

This paper is organized as follows. In Section 2, we briefly review the vCH equation. We propose the cCS scheme for the vCH equation in Section 3. In Sections 4 and 5, we present numerical experiments with various free energies. Finally, conclusions are drawn in Section 6.

# 2. VECTOR-VALUED CAHN–HILLIARD EQUATION

Let c = (c1, . . . , cN )T be the vector-valued concentration field. Clearly the concentration fields satisfy the partition of unity constraint,

c1 + · · · + cN = 1.       (2.1)

Hence admissible states belong to Gibbs N-simplex G :=$\left\{\mathbf{c} \in \mathbb{R}^{N} \mid \sum_{i=1}^{N} c_{i}=1,0 \leq c_{i} \leq 1\right\}$. There are several generalizations of the GL energy functional (1.1) to multi-component systems. In this paper, we adopt the approach used in [5, 10, 12, 18, 19, 20, 21, 22, 23, 27, 31, 32, 33, 34], where the energy functional is defined as follows:

ε(c) :=$\int_{\Omega}\left(\mathcal{F}(\mathbf{c})+\frac{\epsilon^{2}}{2} \sum_{i=1}^{N}\left|\nabla c_{i}\right|^{2}\right)$dx,       (2.2)

referred to as the vGL the energy functional. There are two typical choices for F(c): the polynomial free energy [20, 21, 22, 23, 27, 31, 32, 33, 34]

$\mathcal{F}(\mathbf{c})=\frac{1}{4} \sum_{i=1}^{N} c_{i}^{2}\left(c_{i}-1\right)^{2}$

and logarithmic free energy [5, 10, 12, 18, 19, 32]

$\mathcal{F}(\mathbf{c})=\theta \sum_{i=1}^{N} c_{i} \ln c_{i}+\theta_{c} \sum_{i=1}^{N} \sum_{j=i+1}^{N} c_{i} c_{j}=\theta \sum_{i=1}^{N} c_{i} \ln c_{i}+\frac{\theta_{c}}{2} \sum_{i=1}^{N} c_{i}\left(1-c_{i}\right)$,

where θ and θc are the absolute and critical temperatures, respectively. The vCH equation is a gradient flow for the vGL energy functional (2.2) in the H−1 inner product under the constraint (2.1). In order to ensure (2.1), we use a Lagrange multiplier α(c) [6, 16, 18, 19, 21, 22, 23, 26, 35]. Then, the vCH equation becomes

$\frac{\partial \mathbf{c}}{\partial t}$=∆µ, µ := $\frac{\delta \mathcal{E}}{\delta \mathbf{c}}$= f(c) − ϵ2∆c + α(c)1, x ∈ Ω, 0 < t ≤ T,

ci(x, 0) = $c_{i}^{0}$(x), ∇ci · n = ∇µi · n = 0 on ∂Ω, for i = 1, . . . , N,       (2.3)

where µ = (µ1, . . . , µN )T is the vector-valued chemical potential, µi is the chemical potential of the phase i, $\frac{\delta}{\delta \mathbf{c}}$denotes the variational derivative with respect to c, f(c) = $\frac{\partial \mathcal{F}}{\partial \mathbf{c}}$ = $\left(\frac{\partial \mathcal{F}}{\partial c_{1}}, \ldots, \frac{\partial \mathcal{F}}{\partial c_{N}}\right)^{T}$ = (f(c1), . . . , f(cN ))T , 1 = (1, . . . , 1)T ∈ ℝN , n is a unit normal vector to ∂Ω, and

$\alpha(\mathbf{c})=-\frac{1}{N} \sum_{i=1}^{N} f\left(c_{i}\right)$.

Because the vCH equation (2.3) is of gradient type, the vGL energy functional is nonincreasing in time as the partition of unity constraint holds:

\begin{aligned} \frac{d \mathcal{E}}{d t} &=\int_{\Omega} \sum_{i=1}^{N}\left(\frac{\partial \mathcal{F}}{\partial c_{i}} \frac{\partial c_{i}}{\partial t}+\epsilon^{2} \nabla c_{i} \cdot \nabla \frac{\partial c_{i}}{\partial t}\right) d \mathbf{x}=\int_{\Omega} \sum_{i=1}^{N}\left(f\left(c_{i}\right)-\epsilon^{2} \Delta c_{i}\right) \frac{\partial c_{i}}{\partial t} d \mathbf{x} \\ &=\int_{\Omega} \sum_{i=1}^{N}\left(\mu_{i}-\alpha(\mathbf{c})\right) \frac{\partial c_{i}}{\partial t} d \mathbf{x}=\int_{\Omega} \sum_{i=1}^{N} \mu_{i} \Delta \mu_{i} d \mathbf{x}-\int_{\Omega} \alpha(\mathbf{c}) \frac{\partial}{\partial t} \sum_{i=1}^{N} c_{i} d \mathbf{x} \\ &=-\int_{\Omega} \sum_{i=1}^{N}\left|\nabla \mu_{i}\right|^{2} d \mathbf{x} \leq 0. \end{aligned}

# 3. CONSTRAINED CONVEX SPLITTING SCHEME FOR THE VECTOR-VALUED CAHN–HILLIARD EQUATION

In this section, we propose an unconditionally energy stable scheme for the vCH equation (2.3). The scheme is based on the observation that the vGL energy functional (2.2) can be split into convex and concave parts:

$\mathcal{E}(\mathbf{c})=\mathcal{E}_{c}(\mathbf{c})-\mathcal{E}_{e}(\mathbf{c})=\int_{\Omega}\left(\mathcal{F}_{c}(\mathbf{c})+\frac{\epsilon^{2}}{2} \sum_{i=1}^{N}\left|\nabla c_{i}\right|^{2}\right) d \mathbf{x}-\int_{\Omega} \mathcal{F}_{e}(\mathbf{c})$ dx,

where Fc(c) and −Fe(c) are convex and concave parts of F(c), respectively. Treating Ec(c) implicitly and Ee(c) explicitly, the constrained Convex Splitting (cCS) scheme is obtained:

\begin{aligned} \frac{\mathbf{c}^{n+1}-\mathbf{c}^{n}}{\Delta t} &=\Delta \boldsymbol{\mu}^{n+1} \\ \boldsymbol{\mu}^{n+1} &:=\frac{\delta \mathcal{E}_{c}\left(\mathbf{c}^{n+1}\right)}{\delta \mathbf{c}}-\frac{\delta \mathcal{E}_{e}\left(\mathbf{c}^{n}\right)}{\delta \mathbf{c}} \end{aligned}       (3.1)

$=\mathbf{f}_{c}\left(\mathbf{c}^{n+1}\right)-\epsilon^{2} \Delta \mathbf{c}^{n+1}+\alpha_{c}\left(\mathbf{c}^{n+1}\right) \mathbf{1}-\mathbf{f}_{e}\left(\mathbf{c}^{n}\right)-\alpha_{e}\left(\mathbf{c}^{n}\right) \mathbf{1}$       (3.2)

where fc(c) = $\frac{\partial \mathcal{F}_{c}}{\partial \mathbf{c}}$= (fc(c1), . . . , fc(cN ))T, fe(c) = $\frac{\partial \mathcal{F}_{c}}{\partial \mathbf{c}}$ = = (fe(c1), . . . , fe(cN ))T,

αc(c) = $-\frac{1}{N} \sum_{i=1}^{N} f_{c}\left(c_{i}\right)$, and αe(c) =$-\frac{1}{N} \sum_{i=1}^{N} f_{e}\left(c_{i}\right)$.

Remark 3.1. The choices of αc(cn+1) and αe(cn) in the cCS scheme (3.1) and (3.2) are critical factors to satisfy the partition of unity constraint (2.1) and have the unconditional energy

stability. αc(cn+1) and αe(cn) are obtained by taking the variational derivative of εc(cn+1)

and εc(cn) with respect to c under the constraints $\sum_{i=1}^{N} c_{i}^{n+1}$ = 1 and $\sum_{i=1}^{N} c_{i}^{n}$ = 1, respectively

Lemma 3.2. The cCS scheme is mass conserving.

Proof. Let $c_{i}^{n+1}$(i = 1, . . . , N) be a solution of the cCS scheme. From Eq. (3.1), we have

$\frac{1}{\Delta t} \int_{\Omega}\left(c_{i}^{n+1}-c_{i}^{n}\right)$ dx = $\int_{\Omega} \Delta \mu_{i}^{n+1}$ dx = $\int_{\partial \Omega} \nabla \mu_{i}^{n+1}$· n ds = 0, for i = 1, . . . , N,

where we used the zero Neumann boundary condition for $\mu_{i}^{n+1}$. It follows that $\int_{\Omega} c_{i}^{n+1}$ dx = $\int_{\Omega} c_{i}^{n}$ dx for each i.

Lemma 3.3. The cCS scheme satisfies the constraint (2.1) at any time t n , i.e., $\sum_{i=1}^{N} c_{i}^{n}$ = 1 if an initial condition satisfies $\sum_{i=1}^{N} c_{i}^{0}$ = 1.

Proof. Since $\sum_{i=1}^{N} f_{c}\left(c_{i}^{n+1}\right)$ + Nαc(cn+1) = 0 and $\sum_{i=1}^{N} f_{e}\left(c_{i}^{n}\right)$ + Nαe(cn) = 0, we have from Eqs. (3.1) and (3.2)

$\frac{1}{\Delta t} \sum_{i=1}^{N}\left(c_{i}^{n+1}-c_{i}^{n}\right)=-\epsilon^{2} \Delta^{2} \sum_{i=1}^{N} c_{i}^{n+1}$,i.e., (I + ∆t ϵ22$\sum_{i=1}^{N} c_{i}^{n+1}=\sum_{i=1}^{N} c_{i}^{n}$,       (3.3)

where I denotes the identity operator. Since I + ∆t ϵ22 with a zero Neumann boundary condition for ci is an invertible operator, Eq. (3.3) has a unique solution. Thus, Eq. (3.3)

ensures that $\sum_{i=1}^{N} c_{i}^{n+1}$ = 1 for all n ≥ 0 with the initial condition satisfying $\sum_{i=1}^{N} c_{i}^{0}$ = 1.

Theorem 3.4. The cCS scheme is uniquely solvable for any time step ∆t > 0

Proof. We consider the following functional on $\widetilde{H}$ = {$\mathbf{c} \mid \sum_{i=1}^{N} c_{i}$ = 1, $\int_{\Omega} c_{i}$ dx = $\int_{\Omega} c_{i}^{n}$ dx for i = 1,...,N }:

G(c) = $\frac{1}{2 \Delta t}\left\|\mathbf{c}-\mathbf{c}^{n}\right\|_{H^{-1}}^{2}+\mathcal{E}_{c}(\mathbf{c})-\left(\frac{\delta \mathcal{E}_{e}\left(\mathbf{c}^{n}\right)}{\delta \mathbf{c}}, \mathbf{c}\right)$L2,

where (c, d)L2 =$\int_{\Omega}$ c · d dx = $\int_{\Omega} \sum_{i=1}^{N}$cidi dx.

It may be shown that cn+1$\widetilde{H}$ is the unique minimizer of G(c) if and only if it solves,

for any d ∈ H0$\left\{\mathbf{d} \mid \sum_{i=1}^{N} d_{i}=0, \int_{\Omega} d_{i} d \mathbf{x}=0 \text { for } i=1, \ldots, N\right\}$,

$\left.\frac{d G(\mathbf{c}+s \mathbf{d})}{d s}\right|_{s=0}=\left(\frac{\mathbf{c}-\mathbf{c}^{n}}{\Delta t}, \mathbf{d}\right)_{H^{-1}}+\left(\frac{\delta \mathcal{E}_{c}(\mathbf{c})}{\delta \mathbf{c}}-\frac{\delta \mathcal{E}_{e}\left(\mathbf{c}^{n}\right)}{\delta \mathbf{c}}, \mathbf{d}\right)_{L^{2}}$       (3.4)

$=\left(\frac{\mathrm{c}-\mathrm{c}^{n}}{\Delta t}-\Delta\left(\frac{\delta \mathcal{E}_{c}(\mathrm{c})}{\delta \mathrm{c}}-\frac{\delta \mathcal{E}_{e}\left(\mathrm{c}^{n}\right)}{\delta \mathrm{c}}\right), \mathrm{d}\right)_{H^{-1}}$= 0,       (3.5)

because G(c) is strictly convex by

$\left.\frac{d^{2} G(\mathbf{c}+s \mathbf{d})}{d s^{2}}\right|_{s=0}=\frac{1}{\Delta t}\|\mathbf{d}\|_{H^{-1}}^{2}+\int_{\Omega} \sum_{i=1}^{N}\left(f_{c}^{\prime}\left(c_{i}\right) d_{i}^{2}+\epsilon^{2}\left|\nabla d_{i}\right|^{2}\right)$dx ≥ 0.

Here, the second term on the right-hand side of Eq. (3.4) is obtained as follows:

$\left.\frac{d}{d s} \int_{\Omega}\left(\mathcal{F}_{c}(\mathbf{c}+s \mathbf{d})+\frac{\epsilon^{2}}{2} \sum_{i=1}^{N}\left|\nabla\left(c_{i}+s d_{i}\right)\right|^{2}\right) d \mathbf{x}\right|_{s=0}$

$\int_{\Omega}\left(\mathbf{f}_{c}(\mathbf{c}) \cdot \mathbf{d}+\epsilon^{2} \sum_{i=1}^{N} \nabla c_{i} \cdot \nabla d_{i}\right) d \mathbf{x}=\int_{\Omega}\left(\mathbf{f}_{c}(\mathbf{c})-\epsilon^{2} \Delta \mathbf{c}\right)$· d dx

$\int_{\Omega}\left(\mathbf{f}_{c}(\mathbf{c})-\epsilon^{2} \Delta \mathbf{c}+\alpha_{c}(\mathbf{c}) \mathbf{1}\right) \cdot \mathbf{d} d \mathbf{x}-\int_{\Omega} \alpha_{c}(\mathbf{c})$1 · d dx = $\left(\frac{\delta \mathcal{E}_{c}(\mathbf{c})}{\delta \mathbf{c}}, \mathbf{d}\right)_{L^{2}}$.

And, Eq. (3.5) is true for any d ∈ H0 if and only if the given equation holds

$\frac{\mathbf{c}^{n+1}-\mathbf{c}^{n}}{\Delta t}=\Delta\left(\frac{\delta \mathcal{E}_{c}\left(\mathbf{c}^{n+1}\right)}{\delta \mathbf{c}}-\frac{\delta \mathcal{E}_{e}\left(\mathbf{c}^{n}\right)}{\delta \mathbf{c}}\right)$.

Hence, minimizing the strictly convex functional G(c) is equivalent to solving Eqs. (3.1) and (3.2).

Lemma 3.5. Consider the following convex splitting of the GL energy functional (1.1):

E(c) = Ec(c)−Ee(c) = $\int_{\Omega}\left(F_{c}(c)+\frac{\epsilon^{2}}{2}|\nabla c|^{2}\right)$ dx - $\int_{Ω}$ Fe(c) dx, where $F_{c}^{\prime}$ (c) = fc(c) and

$F_{c}^{\prime}$(c) = fe(c). Then, the convexity of Ec(c) and Ee(c) yields the following inequality:

E(cn+1) − E(cn) ≤ $\int_{\Omega}\left(\frac{\delta E_{c}\left(c^{n+1}\right)}{\delta c}-\frac{\delta E_{e}\left(c^{n}\right)}{\delta c}\right)$(cn+1 − cn ) dx

$\int_{Ω}$ ( fc(cn+1) − ϵ2∆cn+1 − fe(cn)) (cn+1 − cn ) dx.

Proof. We refer to [37].

Theorem 3.6. The cCS scheme with an initial condition satisfying $\sum_{i=1}^{N} c_{i}^{0}$ = 1 is unconditionally energy stable, meaning that for any ∆t > 0,

ε(cn+1) ≤ ε(cn).

Proof. Using Lemma 3.5, we have

ε(cn+1) − ε(cn) = $\sum_{i=1}^{N}\left(E\left(c_{i}^{n+1}\right)-E\left(c_{i}^{n}\right)\right)$

≤ $\sum_{i=1}^{N} \int_{\Omega}\left(f_{c}\left(c_{i}^{n+1}\right)-\epsilon^{2} \Delta c_{i}^{n+1}-f_{e}\left(c_{i}^{n}\right)\right)\left(c_{i}^{n+1}-c_{i}^{n}\right)$dx.

And we obtain from Lemma 3.3

ε(cn+1) − ε(cn) ≤ $\sum_{i=1}^{N} \int_{\Omega}\left(\mu_{i}^{n+1}-\alpha_{c}\left(\mathbf{c}^{n+1}\right)+\alpha_{e}\left(\mathbf{c}^{n}\right)\right)\left(c_{i}^{n+1}-c_{i}^{n}\right)$dx

$\Delta t \sum_{i=1}^{N} \int_{\Omega} \mu_{i}^{n+1} \Delta \mu_{i}^{n+1}$dx

$-\int_{\Omega}\left(\alpha_{c}\left(\mathbf{c}^{n+1}\right)-\alpha_{e}\left(\mathbf{c}^{n}\right)\right) \sum_{i=1}^{N}\left(c_{i}^{n+1}-c_{i}^{n}\right)$dx

$-\Delta t \sum_{i=1}^{N} \int_{\Omega}\left|\nabla \mu_{i}^{n+1}\right|^{2}$dx ≤ 0.

# 4. NUMERICAL EXPERIMENTS WITH THE POLYNOMIAL FREE ENERGY

For numerical tests, we consider the polynomial free energy F(c) = $\frac{1}{4} \sum_{i=1}^{N} c_{i}^{2}$(ci - 1)2 and following splitting:

Fc(c) = $\frac{1}{4} \sum_{i=1}^{N} c_{i}^{2}, \mathcal{F}_{e}(\mathbf{c})=\sum_{i=1}^{N} \Psi\left(c_{i}\right), \Psi(c):=\left\{\begin{array}{ll} 0, & c<0 \\ \frac{1}{4}\left(2 c^{3}-c^{4}\right), & c \in[0,1] \\ \frac{1}{4}(-1+2 c), & c>1 \end{array}\right.$       (4.1)

to have both convexity of Fc(c) and Fe(c) for the extended range of c and easiness of implementation. Then, fc(c) = $1\over2$c and fe(c) = Ψ′(c), and the cCS scheme with Fc(c) and Fe(c) in (4.1) becomes

$\frac{\mathbf{c}^{n+1}-\mathbf{c}^{n}}{\Delta t}=\Delta\left(\frac{1}{2} \mathbf{c}^{n+1}-\epsilon^{2} \Delta \mathbf{c}^{n+1}-\frac{1}{2 N} \sum_{i=1}^{N} c_{i}^{n+1} \mathbf{1}-\mathbf{f}_{e}\left(\mathbf{c}^{n}\right)-\alpha_{e}\left(\mathbf{c}^{n}\right) \mathbf{1}\right)$.       (4.2)

By Lemma 3.3, Eq. (4.2) can be rewritten as follows:

$\frac{\mathbf{c}^{n+1}-\mathbf{c}^{n}}{\Delta t}=\Delta\left(\frac{1}{2} \mathbf{c}^{n+1}-\epsilon^{2} \Delta \mathbf{c}^{n+1}-\mathbf{f}_{e}\left(\mathbf{c}^{n}\right)-\alpha_{e}\left(\mathbf{c}^{n}\right) \mathbf{1}\right)$.

Since fc(c) is linear with respect to c, the scheme allows to solve the vCH equation componentwisely,

$\mathcal{D} c_{i}^{n+1}=b_{i}^{n}$, i.e., $c_{i}^{n+1}=\mathcal{D}^{-1} b_{i}^{n}$, for i = 1, . . . , N,

where D := I − ∆t∆($1\over2$− ϵ2∆) is invertible with a zero Neumann boundary condition for ci and $b_{i}^{n}$ := $c_{i}^{n}$− ∆t∆ ( fe($c_{i}^{n}$)+ αe(cn)) is given explicitly.

We here use the Fourier spectral method for the spatial discretization and the discrete cosine transform in MATLAB is applied for the whole numerical simulations to solve the vCH equation with the zero Neumann boundary condition.

4.1. Convergence test. We demonstrate the convergence of the proposed scheme with the initial conditions

c1(x, 0) = $1\over3$+ 0.01 cos$3\over2$x, c2(x, 0) =$1\over3$+ 0.02 cos x,

c3(x, 0) = 1 − c1(x, 0) − c2(x, 0)

on Ω = [0, 2π]. We set ϵ = 0.25 and compute c(x, t) for 0 < t ≤ 280. The grid size is fixed to h = 2π/128 which provides enough spatial accuracy. In order to estimate the convergence rate with respect to ∆t, simulations are performed by varying ∆t = 2−10 , 2−9 , . . . , 22 . We take the quadruply over-resolved numerical solution as the reference solution. Figures 1 (a) and (b) show the evolution of ε(t) for the reference solution and the relative l2-errors of c(x, 120) (this time is indicated by a dashed line in Fig. 1 (a)) for various time steps, respectively. It is observed that the scheme is first-order accurate in time.

FIGURE 1. (a) Evolution of ε(t) for the reference solution with ϵ = 0:25 and h = 2π/128. (b) Relative l2-errors of c(x, 120) for various time steps.

4.2. Energy stability of the proposed scheme. In order to investigate the energy stability of the proposed scheme, we consider the phase separation of a ternary system with the initial conditions

c1(x, y, 0) = $1\over3$+ rand(x, y), c2(x, y, 0) =$1\over3$+ rand(x, y),

c3(x, y, 0) = 1 − c1(x, y, 0) − c2(x, y, 0)

on Ω = [0, 2π] × [0, 2π]. Here, rand(x, y) is a random number between −0.1 and 0.1, and we use ϵ = 0.1 and h = 2π/128.

Figure 2 shows the evolution of ε(t) using the explicit Euler’s and the proposed scheme with different time steps. The energy curves for the explicit Euler scheme with ∆t = 2−20 and 2−19 are nonincreasing until t = 2−14, whereas the energy curve with ∆t = 2−18 increases rapidly after t = 2−16. As we can see this figure, the explicit Euler’s scheme has a severe time step restriction for energy stability. All the energy curves using the proposed scheme are nonincreasing in time even for sufficiently large time steps. Figure 3 shows the evolution of c(x, y, t) with ∆t = 2−10.

FIGURE 2. Evolution of ε(t) using the explicit Euler’s and the proposed scheme with different time steps.

FIGURE 3. Evolution of c(x, y, t) with ϵ = 0:1, h = 2 π / 128, and Δt = 2-10. In each snapshots, the yellow, green, and blue regions indicate c1, c2, and c3, respectively, and contour lines represent ci = 0:5.

4.3. Efficiency of the proposed scheme. In order to show the efficiency of the proposed scheme, we consider the phase separation of 3, . . . , 10 components (N = 3, . . . , 10). For each N, the initial conditions are chosen as follows: the domain Ω = [0, 2π] × [0, 2π] is partitioned into 40 Voronoi cells and ci is set to 1 on randomly selected Voronoi cells for i = 1, . . . , N. ϵ = 0.05, h = 2π/128, and ∆t = 1/4 are used. Simulations are run until T = 512 and Fig. 4 shows the evolution of c(x, y, t) at t = 0 and 512. Figure 5 presents the CPU time (in seconds, averaged over 10 trials performed on Intel Core i5-7500 CPU at 3.40GHz with 8GB RAM) consumed for N = 3, . . . , 10. The results suggest that the CPU time is almost linear with respect to the number of components N.

FIGURE 4. Evolution of c(x, y, t) at t = 0 (top) and 512 (bottom) with ϵ = 0:05, h = 2π / 128, and Δt = 1=4. Columns 1–4 correspond to N = 3, 5, 8, and 10, respectively. In the top, the 40 Voronoi cells are represented by red dotted lines.

FIGURE 5. CPU time versus the number of components. Each simulation is run until T = 512. Each line segment is obtained by least squares ﬁtting of all points.

4.4. Phase separation of a four-component mixture in 3D. We solve the vCH equation on Ω = [0, 2π] × [0, 2π] × [0, 2π] with ϵ = 0.1, h = 2π/64, and ∆t = 1/4. The initial conditions are

ci(x, y, z, 0) = $1\over4$ + i · rand(x, y, z), for i = 1, 2, 3,

c4(x, y, z, 0) = 1 - $\sum_{i=1}^{3} c_{i}$(x, y, z, 0),

where rand(x, y, z) is a random number between −0.01 and 0.01 at the grid points. Figures 6 and 7 show the evolution of c(x, y, z, t) and its energy, respectively. We observe that the randomly perturbed constant concentration fields evolve to many small structures and then to single-interface structures as the energy is dissipated in time.

FIGURE 6. Evolution of c(x, y, z, t) with ϵ = 0:1, h = 2π / 64, and Δt = 1=4. In each snapshots, the red, green, blue, and yellow regions indicate c1, c2, c3, and c4, respectively.

FIGURE 7. Evolution of ϵ(t).

# 5. NUMERICAL EXPERIMENTS WITH THE LOGARITHMIC FREE ENERGY In this section, we consider the logarithmic free energy

F(c) = θ$\sum_{i=1}^{N} c_{i}$ln ci +$\theta_{c} \sum_{i=1}^{N} \sum_{j=i+1}^{N} c_{i} c_{j}$ = $\theta \sum_{i=1}^{N} c_{i}$ln ci$θc\over2$$\sum_{i=1}^{N} c_{i}$(1 − ci).

Unlike the polynomial free energy, a nonlinear convex splitting is a natural choice for the logarithmic free energy:

Fc(c) = θ $\sum_{i=1}^{N} c_{i}$ ln ci , Fe(c) = $-\frac{\theta_{c}}{2} \sum_{i=1}^{N} c_i$(1 − ci).

Then, fc(c) = θ ln c and fe(c) = θcc. In the case of the logarithmic free energy, there is a numerical difficulty associated with the singularity as each ci approaches zero. In order to avoid this, we apply a regularization to the logarithmic function, i.e., for a small positive number δ, we define

lnδ c := $\left\{\begin{array}{ll} \ln c, & \text { if } c \geq \delta, \\ p(c)=-\frac{1}{2 \delta^{2}} c^{2}+\frac{2}{\delta} c+\ln \delta-\frac{3}{2}, & \text { otherwise, } \end{array}\right.$

where the quadratic polynomial p(c) matches the values of zeroth, first, and second derivatives of the logarithmic function at c = δ [3, 4].

The nonlinearity of the scheme comes from fc($c_{i}^{n+1}$) and αc(cn+1) and these can be handled using a Newton-type linearization [24, 28]

$f_{c}\left(c_{i}^{n, m+1}\right) \approx f_{c}\left(c_{i}^{n, m}\right)+f_{c}^{\prime}\left(c_{i}^{n, m}\right)\left(c_{i}^{n, m+1}-c_{i}^{n, m}\right)$

$\alpha_{c}\left(\mathbf{c}^{n, m+1}\right) \approx \alpha_{c}\left(\mathbf{c}^{n, m}\right)+\frac{\partial \alpha_{c}\left(\mathbf{c}^{n, m}\right)^{T}}{\partial \mathbf{c}}\left(\mathbf{c}^{n, m+1}-\mathbf{c}^{n, m}\right)$

$=-\frac{1}{N} \sum_{i=1}^{N}\left(f_{c}\left(c_{i}^{n, m}\right)+f_{c}^{\prime}\left(c_{i}^{n, m}\right)\left(c_{i}^{n, m+1}-c_{i}^{n, m}\right)\right)$

for m = 0, 1, . . .. We then develop a Newton-type fixed point iteration method for the cCS scheme as

$\left(\begin{array}{cccc} \mathcal{D}_{1}+\mathcal{A}_{1} & \mathcal{A}_{2} & \cdots & \mathcal{A}_{N} \\ \mathcal{A}_{1} & \mathcal{D}_{2}+\mathcal{A}_{2} & \cdots & \mathcal{A}_{N} \\ \vdots & \vdots & \ddots & \vdots \\ \mathcal{A}_{1} & \mathcal{A}_{2} & \cdots & \mathcal{D}_{N}+\mathcal{A}_{N} \end{array}\right)\left(\begin{array}{c} c_{1}^{n, m+1}-c_{1}^{n, m} \\ c_{2}^{n, m+1}-c_{2}^{n, m} \\ \vdots \\ c_{N}^{n, m+1}-c_{N}^{n, m} \end{array}\right)=\left(\begin{array}{c} b_{1}^{n, m} \\ b_{2}^{n, m} \\ \vdots \\ b_{N}^{n, m} \end{array}\right)$  ,     (5.1)

where cn,0 = cn,

Di = I − ∆t ∆($f_{c}^{\prime}\left(c_{i}^{n, m}\right)$− ϵ 2∆ ) , Ai = −∆t ∆ $\left(-\frac{1}{N} f_{c}^{\prime}\left(c_{i}^{n, m}\right)\right)$,

$b_{i}^{n, m}=c_{i}^{n}-c_{i}^{n, m}+\Delta t \Delta\left(f_{c}\left(c_{i}^{n, m}\right)-\epsilon^{2} \Delta c_{i}^{n, m}+\alpha_{c}\left(\mathbf{c}^{n, m}\right)-f_{e}\left(c_{i}^{n}\right)-\alpha_{e}\left(\mathbf{c}^{n}\right)\right)$,

for i = 1, . . . , N, and we set

cn+1 = cn,m+1

if a relative l2-norm of the consecutive error$\frac{\left\|\mathbf{c}^{n, m+1}-\mathbf{c}^{n, m}\right\|_{2}}{\left\|\mathbf{c}^{n, m}\right\|_{2}}$is less than a tolerance tol. In this paper, the biconjugate gradient (BICG) method is used to solve the system (5.1) and we use the following preconditioner P to accelerate the convergence speed of the BICG algorithm:

$P=\left(\begin{array}{cccc} \overline{\mathcal{D}}_{1} & 0 & \cdots & 0 \\ 0 & \overline{\mathcal{D}}_{2} & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & \overline{\mathcal{D}}_{N} \end{array}\right)$,

where  $\overline{\mathcal{D}}_{i}$ = I − ∆t ∆$\left(\overline{f_{c}^{\prime}\left(c_{i}^{n, m}\right)}-\epsilon^{2} \Delta\right)$and $\overline{f_{c}^{\prime}\left(c_{i}^{n, m}\right)}$is the average value of $f_{c}^{\prime}\left(c_{i}^{n, m}\right)$.The stopping criterion for the BICG iteration is that the relative residual norm is less than tol.

5.1. Robustness of the nonlinear solver and convergence test. In order to show the robustness of the nonlinear solver and the necessity of the preconditioner, we count the number of nonlinear and BICG iterations with the initial conditions

$c_{1}(x, 0)=\frac{1}{3}+0.01 \cos \frac{3}{2} x, \quad c_{2}(x, 0)=\frac{1}{3}+0.02 \cos x$,

$c_{3}(x, 0)=1-c_{1}(x, 0)-c_{2}(x, 0)$

on Ω = [0, 2π], θ = 0.3, θc = 1, ϵ = 0.25, h = 2π/128, and tol = 10−10. Figure 8 shows the number of BICG iterations without and with the preconditioner during the simulation time 0 < t = n∆t ≤ 280 for different time steps. As shown in Fig. 8, the BICG iterations were remarkably reduced (about 100 times) by using the preconditioner. Figures 9 (a) and (b)

FIGURE 8. Number of BICG iterations without and with the preconditioner for different time steps.

FIGURE 9. Number of (a) nonlinear and (b) BICG iterations (with the precon-ditioner) averaged over the simulation time.

show the number of nonlinear and BICG iterations (with the preconditioner) averaged over the simulation time, respectively. 2–4 nonlinear iterations (on average) were involved in proceeding to the next time level. We believe that such a fast iterative convergence can be achieved since the successive iteration (5.1) is a Newton-type fixed point iteration method. And 6–42 BICG iterations (on average) were involved in proceeding to the next time level owing to the preconditioner.

Next, we vary ∆t = 2−10, 2−9, . . . , 22 to estimate the convergence rate with respect to ∆t for the logarithmic free energy. We take the quadruply over-resolved numerical solution as the reference solution. Figures 10 (a) and (b) show the evolution of ε(t) for the reference solution and the relative l2-errors of c(x, 120) (this time is indicated by a dashed line in Fig. 10 (a)) for various time steps, respectively. It is observed that the scheme is also first-order accurate in time for the logarithmic free energy.

FIGURE 10. (a) Evolution of ϵ(t) for the reference solution with θ = 0:3, θc = 1, ϵ = 0:25, h = 2π / 128. (b) Relative l2-errors of c(x, 120) for various time steps.

5.2. Energy stability of the proposed scheme. In order to investigate the energy stability for the logarithmic free energy, we take the initial conditions as

c1(x, y, 0) = $1\over3$+ rand(x, y), c2(x, y, 0) = $1\over3$+ rand(x, y),

c3(x, y, 0) = 1 − c1(x, y, 0) − c2(x, y, 0)

on Ω = [0, 2π] × [0, 2π]. Here, rand(x, y) is a random number between −0.1 and 0.1, and we use θ = 0.3, θc = 1, ϵ = 0.1, h = 2π/128, and tol = 10−6 . Figure 11 shows the evolution of ε(t) with different time steps. All the energy curves are also nonincreasing in time for the logarithmic free energy. Figure 12 shows the evolution of c(x, y, t) with ∆t = 2−2 .

FIGURE 11. Evolution of ϵ(t) with different time steps.

FIGURE 12. Evolution of c(x, y, t) with θ = 0:3, θc = 1, ϵ = 0:1, h = 2π / 128, and Δt = 2-2. In each snapshots, the yellow, green, and blue regions indicate c1, c2, and c3, respectively, and contour lines represent ci = 0:45, 0.5, and 0.55.

# 6. CONCLUSIONS

In this paper, we proposed the cCS scheme for the vCH equation and proved its unconditional energy stability. For the polynomial free energy and linear convex splitting, we confirmed that the scheme is first-order accurate in time and unconditionally energy stable. Owing to the linear convex splitting, we solved the vCH equation efficiently (the CPU time was almost linear with respect to the number of components N). For the logarithmic free energy and nonlinear convex splitting, we showed the robustness of the nonlinear solver and the necessity of the preconditioner. And we also demonstrated that the scheme is first-order accurate in time and unconditionally energy stable.

We note that order of time accuracy of the cCS scheme can be improved by various approaches. One of them is to combine with an s-stage implicit–explicit Runge–Kutta method [1, 29, 30] and extension of the cCS scheme to high-order time accuracy can be considered as the scope of future study.

#### Acknowledgement

Supported by : National Research Foundation (NRF)

#### References

1. U.M. Ascher, S.J. Ruuth, and R.J. Spiteri, Implicit-explicit Runge-Kutta methods for time-dependent partial differential equations, Appl. Numer. Math. 25 (1997), 151-167. https://doi.org/10.1016/S0168-9274(97)00056-1
2. V.E. Badalassi, H.D. Ceniceros, and S. Banerjee, Computation of multiphase systems with phase field models, J. Comput. Phys. 190 (2003), 371-397. https://doi.org/10.1016/S0021-9991(03)00280-8
3. J.W. Barrett and J.F. Blowey, An error bound for the finite element approximation of the Cahn-Hilliard equation with logarithmic free energy, Numer. Math. 72 (1995), 1-20. https://doi.org/10.1007/s002110050157
4. J.W. Barrett and J.F. Blowey, An error bound for the finite element approximation of a model for phase separation of a multi-component alloy, IMA J. Numer. Anal. 16 (1996), 257-287. https://doi.org/10.1093/imanum/16.2.257
5. J.F. Blowey, M.I.M. Copetti, and C.M. Elliott, Numerical analysis of a model for phase separation of a multicomponent alloy, IMA J. Numer. Anal. 16 (1996), 111-139. https://doi.org/10.1093/imanum/16.1.111
6. F. Boyer and C. Lapuerta, Study of a three component Cahn-Hilliard flow model, ESAIM: M2AN 40 (2006), 653-687. https://doi.org/10.1051/m2an:2006028
7. F. Boyer and S. Minjeaud, Hierarchy of consistent n-component Cahn-Hilliard systems, Math. Models Meth. Appl. Sci. 24 (2014), 2885-2928. https://doi.org/10.1142/S0218202514500407
8. J.W. Cahn and J.E. Hilliard, Free energy of a nonuniform system. I. Interfacial free energy, J. Chem. Phys. 28 (1958), 258-267. https://doi.org/10.1063/1.1744102
9. L.-Q. Chen, Phase-field models for microstructure evolution, Annu. Rev. Mater. Res. 32 (2000), 113-140. https://doi.org/10.1146/annurev.matsci.32.112001.132041
10. J.P. Desi, H.H. Edrees, J.J. Price, E. Sander, and T. Wanner, The dynamics of nucleation in stochastic Cahn-Morral systems, SIAM J. Appl. Dyn. Syst. 10 (2011), 707-743. https://doi.org/10.1137/100801378
11. S. Dong, An efficient algorithm for incompressible N-phase flows, J. Comput. Phys. 276 (2014), 691-728. https://doi.org/10.1016/j.jcp.2014.08.002
12. C.M. Elliott and S. Luckhaus, A generalised diffusion equation for phase separation of a multi-component mixture with interfacial free energy, IMA Preprint Series 887, 1991.
13. C.M. Elliott and A.M. Stuart, The global dynamics of discrete semilinear parabolic equations, SIAM J. Numer. Anal. 30 (1993), 1622-1663. https://doi.org/10.1137/0730084
14. D.J. Eyre, Unconditionally gradient stable time marching the Cahn-Hilliard equation, MRS Proc. 529 (1998), 39-46. https://doi.org/10.1557/PROC-529-39
15. D. de Fontaine, A computer simulation of the evolution of coherent composition variations in solid solutions, Ph.D. Thesis, Northwestern University, 1967.
16. H. Garcke, B. Nestler, and B. Stoth, On anisotropic order parameter models for multi-phase systems and their sharp interface limits, Phys. D 115 (1998), 87-108. https://doi.org/10.1016/S0167-2789(97)00227-3
17. Z. Guan, C. Wang, and S.M. Wise, A convergent convex splitting scheme for the periodic nonlocal Cahn- Hilliard equation, Numer. Math. 128 (2014), 377-406. https://doi.org/10.1007/s00211-014-0608-2
18. D. Jeong and J. Kim, A practical numerical scheme for the ternary Cahn-Hilliard system with a logarithm free energy, Phys. A 442 (2016), 510-522. https://doi.org/10.1016/j.physa.2015.09.038
19. D. Jeong and J. Kim, Practical estimation of a splitting parameter for a spectral method for the ternary Cahn-Hilliard system with a logarithmic free energy, Math. Meth. Appl. Sci. 40 (2017), 1734-1745. https://doi.org/10.1002/mma.4093
20. J. Kim, K. Kang, and J. Lowengrub, Conservative multigrid methods for ternary Cahn-Hilliard systems, Comm. Math. Sci. 2 (2004), 53-77. https://doi.org/10.4310/CMS.2004.v2.n1.a4
21. H.G. Lee, J.-W. Choi, and J. Kim, A practically unconditionally gradient stable scheme for the N-component Cahn-Hilliard system, Phys. A 391 (2012), 1009-1019. https://doi.org/10.1016/j.physa.2011.11.032
22. H.G. Lee and J. Kim, A second-order accurate non-linear difference scheme for the N-component Cahn- Hilliard system, Phys. A 387 (2008), 4787-4799. https://doi.org/10.1016/j.physa.2008.03.023
23. H.G. Lee and J. Kim, An efficient and accurate numerical algorithm for the vector-valued Allen-Cahn equations, Comput. Phys. Commun. 183 (2012), 2107-2115. https://doi.org/10.1016/j.cpc.2012.05.013
24. H.G. Lee, J. Shin, and J.-Y. Lee, First- and second-order energy stable methods for the modified phase field crystal equation, Comput. Methods Appl. Mech. Engrg. 321 (2017), 1-17. https://doi.org/10.1016/j.cma.2017.03.033
25. J.E. Morral and J.W. Cahn, Spinodal decomposition in ternary systems, Acta Metall. 19 (1971), 1037-1045. https://doi.org/10.1016/0001-6160(71)90036-8
26. B. Nestler and A.A. Wheeler, A multi-phase-field model of eutectic and peritectic alloys: numerical simulation of growth structures, Phys. D 138 (2000), 114-133. https://doi.org/10.1016/S0167-2789(99)00184-0
27. E. Oudet, Approximation of partitions of least perimeter by $\Gamma$-convergence: around Kelvin's conjecture, Exp. Math. 20 (2011), 260-270. https://doi.org/10.1080/10586458.2011.565233
28. J. Shin, H.G. Lee, and J.-Y. Lee, First and second order numerical methods based on a new convex splitting for phase-field crystal equation, J. Comput. Phys. 327 (2016), 519-542. https://doi.org/10.1016/j.jcp.2016.09.053
29. J. Shin, H.G. Lee, and J.-Y. Lee, Convex Splitting Runge-Kutta methods for phase-field models, Comput. Math. Appl. 73 (2017), 2388-2403. https://doi.org/10.1016/j.camwa.2017.04.004
30. J. Shin, H.G. Lee, and J.-Y. Lee, Unconditionally stable methods for gradient flow using Convex Splitting Runge-Kutta scheme, J. Comput. Phys. 347 (2017), 367-381. https://doi.org/10.1016/j.jcp.2017.07.006
31. R. Tavakoli, Computationally efficient approach for the minimization of volume constrained vector-valued Ginzburg-Landau energy functional, J. Comput. Phys. 295 (2015), 355-378. https://doi.org/10.1016/j.jcp.2015.04.020
32. R. Tavakoli, Unconditionally energy stable time stepping scheme for Cahn-Morral equation: Application to multi-component spinodal decomposition and optimal space tiling, J. Comput. Phys. 304 (2016), 441-464. https://doi.org/10.1016/j.jcp.2015.10.018
33. G.I. Toth, T. Pusztai, and L. Granasy, Consistent multiphase-field theory for interface driven multidomain dynamics, Phys. Rev. B 92 (2015), 184105. https://doi.org/10.1103/PhysRevB.92.184105
34. G.I. Toth, M. Zarifi, and B. Kvamme, Phase-field theory of multicomponent incompressible Cahn-Hilliard liquids, Phys. Rev. E 93 (2016), 013126. https://doi.org/10.1103/PhysRevE.93.013126
35. L. Vanherpe, F. Wendler, B. Nestler, and S. Vandewalle, A multigrid solver for phase field simulation of microstructure evolution, Math. Comput. Simul. 80 (2010), 1438-1448. https://doi.org/10.1016/j.matcom.2009.10.007
36. C.Wang, X.Wang, and S.M.Wise, Unconditionally stable schemes for equations of thin film epitaxy, Discrete Cont. Dyn. S. 28 (2010), 405-423. https://doi.org/10.3934/dcds.2010.28.405
37. C.Wang and S.M.Wise, An energy stable and convergent finite-difference scheme for the modified phase field crystal equation, SIAM J. Numer. Anal. 49 (2011), 945-969. https://doi.org/10.1137/090752675
38. S.M. Wise, C. Wang, and J.S. Lowengrub, An energy-stable and convergent finite-difference scheme for the phase field crystal equation, SIAM J. Numer. Anal. 47 (2009), 2269-2288. https://doi.org/10.1137/080738143