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

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⁻¹ 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) = E_c(c) − E_e(c) = \(\int_{\Omega}\left(F_{c}(c)+\frac{\epsilon^{2}}{2}|\nabla c|^{2}\right) d \mathbf{x}-\int_{\Omega}\) F_e(c) dx,

where F_c(c) and −F_e(c) are convex and concave parts of F(c), respectively. And E_c(c) and E_e(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 c₁, , c_N−1, where c_i is the concentration field of the phase i and N is the number of phases, and enforced cN = 1−\(\sum_{i=1}^{N-1}\)c_i to satisfy the constraint at the next time level. The author in [32] used the Schur complement method to solve multi-component systems for c₁, . . . , c_N−1 and the constraint c₁ + · · · + c_N = 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 = (c₁, . . . , c_N )^T be the vector-valued concentration field. Clearly the concentration fields satisfy the partition of unity constraint,

c₁ + · · · + c_N = 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⁻¹ 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) − ϵ²∆c + α(c)1, x ∈ Ω, 0 < t ≤ T, 

c_i(x, 0) = \(c_{i}^{0}\)(x), ∇c_i · n = ∇µi · n = 0 on ∂Ω, for i = 1, . . . , N,       (2.3)

where µ = (µ₁, . . . , µ_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(c₁), . . . , f(c_N ))^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 F_c(c) and −F_e(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 f_c(c) = \(\frac{\partial \mathcal{F}_{c}}{\partial \mathbf{c}}\)= (f_c(c₁), . . . , f_c(c_N ))^T, f_e(c) = \(\frac{\partial \mathcal{F}_{c}}{\partial \mathbf{c}}\) = = (f_e(c₁), . . . , f_e(c_N ))^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(cⁿ⁺¹) and α_e(cⁿ) 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(cⁿ⁺¹) and α_e(cⁿ) are obtained by taking the variational derivative of ε_c(cⁿ⁺¹)

and ε_c(cⁿ) 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(cⁿ⁺¹) = 0 and \(\sum_{i=1}^{N} f_{e}\left(c_{i}^{n}\right)\) + Nα_e(cⁿ) = 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 ϵ²∆² ) \(\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 ϵ²∆² with a zero Neumann boundary condition for c_i 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}\)c_id_i dx.

It may be shown that cⁿ⁺¹ ∈ \(\widetilde{H}\) is the unique minimizer of G(c) if and only if it solves,

for any d ∈ H₀ = \(\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 ∈ H₀ 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) = E_c(c)−E_e(c) = \(\int_{\Omega}\left(F_{c}(c)+\frac{\epsilon^{2}}{2}|\nabla c|^{2}\right)\) dx - \(\int_{Ω}\) F_e(c) dx, where \(F_{c}^{\prime}\) (c) = f_c(c) and

\(F_{c}^{\prime}\)(c) = f_e(c). Then, the convexity of E_c(c) and E_e(c) yields the following inequality:

E(cⁿ⁺¹) − E(cⁿ) ≤ \(\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)\)(cⁿ⁺¹ − cⁿ ) dx

= \(\int_{Ω}\) ( f_c(cⁿ⁺¹) − ϵ²∆cⁿ⁺¹ − f_e(cⁿ)) (cⁿ⁺¹ − cⁿ ) 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,

ε(cⁿ⁺¹) ≤ ε(cⁿ).

Proof. Using Lemma 3.5, we have

ε(cⁿ⁺¹) − ε(cⁿ) = \(\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

ε(cⁿ⁺¹) − ε(cⁿ) ≤ \(\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}\)(c_i - 1)² and following splitting:

F_c(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 F_c(c) and F_e(c) for the extended range of c and easiness of implementation. Then, f_c(c) = \(1\over2\)c and f_e(c) = Ψ′(c), and the cCS scheme with F_c(c) and F_e(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 f_c(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\)− ϵ²∆) is invertible with a zero Neumann boundary condition for c_i and \(b_{i}^{n}\) := \(c_{i}^{n}\)− ∆t∆ ( fe(\(c_{i}^{n}\))+ α_e(cⁿ)) 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

c₁(x, 0) = \(1\over3\)+ 0.01 cos\(3\over2\)x, c₂(x, 0) =\(1\over3\)+ 0.02 cos x,

c₃(x, 0) = 1 − c₁(x, 0) − c₂(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⁻¹⁰ , 2⁻⁹ , . . . , 2² . 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 l₂-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 l₂-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

c₁(x, y, 0) = \(1\over3\)+ rand(x, y), c₂(x, y, 0) =\(1\over3\)+ rand(x, y),

c₃(x, y, 0) = 1 − c₁(x, y, 0) − c₂(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⁻²⁰ and 2⁻¹⁹ are nonincreasing until t = 2⁻¹⁴, whereas the energy curve with ∆t = 2⁻¹⁸ increases rapidly after t = 2⁻¹⁶. 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⁻¹⁰.

 

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 c₁, c₂, and c₃, respectively, and contour lines represent c_i = 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 c_i 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 fitting 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

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

c₄(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 c₁, c₂, c₃, and c₄, 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 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 + \(θc\over2\)\(\sum_{i=1}^{N} c_{i}\)(1 − c_i).

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

F_c(c) = θ \(\sum_{i=1}^{N} c_{i}\) ln c_i , F_e(c) = \(-\frac{\theta_{c}}{2} \sum_{i=1}^{N} c_i\)(1 − c_i).

Then, f_c(c) = θ ln c and f_e(c) = θcc. In the case of the logarithmic free energy, there is a numerical difficulty associated with the singularity as each c_i 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 f_c(\(c_{i}^{n+1}\)) and α_c(cⁿ⁺¹) 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 c^n,0 = cⁿ,

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

cⁿ⁺¹ = c^n,m+1

if a relative l₂-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⁻¹⁰. 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⁻¹⁰, 2⁻⁹, . . . , 2² 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 l₂-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 l₂-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

c₁(x, y, 0) = \(1\over3\)+ rand(x, y), c₂(x, y, 0) = \(1\over3\)+ rand(x, y),

c₃(x, y, 0) = 1 − c₁(x, y, 0) − c₂(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⁻⁶ . 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⁻² .

 

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 c₁, c₂, and c₃, respectively, and contour lines represent c_i = 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.