A mass conservative and energy stable scheme for the conservative Allen–Cahn type Ohta–Kawasaki model for diblock copolymers

1.   Introduction

The Ohta–Kawasaki (OK) model was introduced to describe phase separation in diblock copolymers ^[1], and has been applied to various systems, including condensed matter and biological systems ^[2,3]. The OK equation,

is the $H^{-1}$-gradient flow of the following energy functional ^[1,4]:

where $\phi$ is the order parameter, $F(\phi) = \frac{1}{4} (\phi^2-1)^2$ is the double-well potential, $\epsilon, \alpha > 0$ are constants, $M > 0$ is a mobility, $\mu$ is the chemical potential, $\frac{\delta}{\delta \phi}$ denotes the variational derivative, $f(\phi) = F'(\phi)$, and the boundary conditions for $\phi$ and $\mu$ are considered periodic in all spatial directions. Additionally, $\psi$ is given by the following solution of the periodic boundary value problem:

where $\bar \phi : = \frac{1}{|\Omega|} \int_\Omega \phi \: d{\mathbf{x}}$.

Recently, the conservative Allen–Cahn type Ohta–Kawasaki (CAC-OK) equation was introduced to reformulate the OK equation ^[5]:

where $\frac{1}{|\Omega|} \int_\Omega f(\phi) \: d{\mathbf{x}}$ is the nonlocal Lagrange multiplier ^[6] (see also ^[7] and its references for theoretical contributions and modeling issues for the conservative Allen–Cahn type model). The CAC-OK equation has the following mass conservation and energy decay properties:

and

Although the CAC-OK equation is of a lower-order than the OK equation, there is one problem: how to discretize $f(\phi)$ and $\frac{1}{|\Omega|} \int_\Omega f(\phi) \: d{\mathbf{x}}$ in time to preserve the mass conservation and energy decay properties without losing the efficiency and accuracy. In ^[5], the scalar auxiliary variable, $u(t) = \sqrt{ \int_\Omega F(\phi) \: d{\mathbf{x}} + C}$, was introduced to redefine the energy functional and reformulate the CAC-OK equation, where $C$ is a constant such that $\int_\Omega F(\phi) \: d{\mathbf{x}} + C > 0$. Moreover, the second-order backward differentiation formula was used to discretize the reformulated system, and the second-order stabilization term, $S(\phi^{n+1} - 2\phi^n + \phi^{n-1})$, was added to improve the energy stability. The aim of this paper is to present a linear, second-order, mass conservative, and energy stable scheme based on the Crank–Nicolson formula for the CAC-OK equation. Here, we shall restrict our attention to $F(\phi)$ that satisfies the following condition: there exists a constant $L > 0$ such that

To satisfy (1.4), we adopt the following truncated double-well potential:

where $p > 0$ is a constant ^[8]. In the scheme, $f(\phi)$ and $\frac{1}{|\Omega|} \int_\Omega f(\phi) \: d{\mathbf{x}}$ are explicitly treated, which make the scheme linear, and the energy stability is guaranteed by adding two second-order stabilization terms. We prove that the scheme is mass conservative and energy stable. Moreover, the scheme can be easily implemented within a few lines of MATLAB code.

The remainder of this paper is organized as follows: we design the numerical scheme and show its mass conservation and energy stability analytically in Section 2 and numerically in Section 3; conclusions are given in Section 4; and the MATLAB code for the numerical scheme is given in the Appendix.

2.   Mass conservative and energy stable scheme

We design the numerical scheme for the CAC-OK equation based on the Crank–Nicolson formula:

where $\phi^{*, n+ \frac{1}{2}} = \frac{3\phi^n-\phi^{n-1}}{2}$, $\phi^{n+ \frac{1}{2}} = \frac{\phi^{n+1} + \phi^n}{2}$, $\psi^{n+ \frac{1}{2}} = \frac{\psi^{n+1} + \psi^n}{2}$, $\delta_t \phi^{n+1} = \phi^{n+1} - \phi^n$, $A, B > 0$ are stabilization parameters, and $\phi^{-1} \equiv \phi^0$.

Theorem 1. The scheme (2.1) is mass conserving.

Proof. Suppose that the scheme (2.1) has a solution. From Equation (2.1), we have the following:

where $(\cdot, \cdot)$ denotes the $L^2$-inner product on $\Omega$. This gives the following relation:

With $\phi^{-1} \equiv \phi^0$, i.e., $\left(\delta_t \phi^0, {\mathbf{1}} \right) = 0$, the relation ensures that

for all $n \geq 0$. □

Theorem 2. The scheme (2.1) with $A \geq \frac{ML^2}{16}$ and $B \geq \frac{L}{2}$ satisfies the following discrete energy decay law:

where

Proof. By simple calculations, we have

and

To handle $f (\phi^{*, n+ \frac{1}{2}})$, we expand $F(\phi^{n+1})$ and $F(\phi^n)$ at $\phi^{*, n+ \frac{1}{2}}$ as follows:

where $\xi^{n+1}$ is a function that is pointwise bounded between $\phi^{n+1}$ and $\phi^{*, n+ \frac{1}{2}}$, and $\xi^n$ is a function that is pointwise bounded between $\phi^n$ and $\phi^{*, n+ \frac{1}{2}}$. Then, we obtain the following:

where we used $\max_{\phi \in \mathbb R} |f'(\phi)| \leq L$. Using (2.3)–(2.6), we have the following:

The right-hand side of the above inequality is less than or equal to zero under $A \geq \frac{ML^2}{16}$ and $B \geq \frac{L}{2}$, and it follows that $\tilde {\mathcal E}^{n+1} \leq \tilde {\mathcal E}^n$. □

Remark 1. The scheme (2.1) is one of the multi-step schemes and its energy stability relies on (1.4). Multi-stage schemes, such as Runge–Kutta, without (1.4) that are applicable to the CAC-OK equation can be found in Refs. ^[9,10,11].

3.   Numerical experiments

3.1.   Numerical implementation

The scheme (2.1) can be rewritten as follows:

For this equation with the periodic boundary condition, we can use the Fourier spectral method ^{[12,13,14,15,16]} to achieve efficient computations. Numerical experiments are implemented in MATLAB, and the MATLAB functions $\texttt{fft2}$, $\texttt{ifft2}$, $\texttt{fftn}$, and $\texttt{ifftn}$ of cost $O(N^d \log N^d)$ are used for the Fourier and inverse Fourier transforms, where $N$ is the number of subintervals in one spatial dimension and $d$ is the number of spatial dimensions.

Table 1 provides the grid size for each numerical experiment.

Unless otherwise stated, we set $p = 1$, $L = 3p^2-1$, $A = \frac{ML^2}{16}$, and $B = \frac{L}{2}$.

3.2.   Efficiency, accuracy, mass conservation, and energy stability tests

First, we check the efficiency and accuracy of the proposed scheme with an initial condition as follows:

where $(x_1, y_1, r_1) = (\pi-0.8, \pi, 1.4)$ and $(x_2, y_2, r_2) = (\pi+1.7, \pi, 0.5)$. We consider the low ($M = 1$, $\alpha = 0.01$) and high ($M = 100$, $\alpha = 100$) stiffness cases, set $\epsilon = 0.06$, and evolve $\phi(x, y, t)$ for $0 < t \leq 10$.

The CPU times (performed on Intel Core i5-7500 CPU at 3.40GHz with 8GB RAM) consumed using the scheme with $\Delta t = \frac{2^{-8}}{M}, \frac{2^{-7}}{M}, \ldots, \frac{1}{M}$ are shown in Figure 1 (a). Additionally, the relative $l_2$-errors of $\phi(x, y, 10)$ are shown in Figure 1 (b), where the error is calculated by comparing with the reference solution using $\Delta t = \frac{2^{-10}}{M}$. It is observed that the CPU time is almost linear with respect to the number of steps, and the convergence order of the scheme is 2.

To see the computational cost with respect to the spatial resolution of the scheme, we fix the time step to $\Delta t = \frac{2^{-4}}{M}$ and vary $N = 128,192,256,384,512$, and $768$. Figure 2 shows the CPU times (performed on Intel Core i5-14400 CPU at 2.50GHz with 8GB RAM) consumed using the scheme. The results indicate that the computational cost is roughly $O(N^2 \log N^2)$.

To explore the influence of the stabilization terms on the accuracy of the scheme, we take various $A$ and $B$ values for the low stiffness case. Figure 3 shows the relative $l_2$-errors of $\phi(x, y, 10)$ with different $A$ and $B$ values. The results indicate the following: (ⅰ) the convergence order of the scheme is still 2 regardless of $A$ and $B$; and (ⅱ) the convergence constant is affected by $A$ and $B$.

Next, we verify the mass conservation and energy stability of the scheme with larger time steps. Figures 4 (a) and (b) show the evolution of $\int_\Omega \left(\phi(x, y, t)- \phi(x, y, 0) \right) dxdy$ and $\tilde {\mathcal E}$ defined in (2.2) with $\Delta t = 2^{-3}, 2^{-2}, \ldots, 2^3$, respectively. As proved by Theorems 1 and 2, the masses are conserving and the energies do not increase in time.

3.3.   Effect of average concentration and total chain length

To investigate the effect of the average concentration and total chain length, we use an initial condition as follows:

where $rand(x, y)$ is a random number between $-0.001$ and $0.001$ at the grid points. We set $M = 1$, $\epsilon = 0.02$, and $\Delta t = \frac{1}{8}$, and evolve $\phi(x, y, t)$ for $0 < t \leq 256$.

With $\alpha = 300000$, Figure 5 shows the evolution of $\phi(x, y, t)$ for $\bar \phi = 0$, 0.1, 0.2, and 0.3. When $\bar \phi = 0$, the final solution forms the stripe pattern. We see the coexistence of short discontinuous stripe and hexagonal spot patterns when $\bar \phi = 0.1$ and 0.2. When $\bar \phi = 0.3$, the hexagonal spot pattern is dominant. These results are consistent with those in ^[17].

Next, Figure 6 shows $\phi(x, y, 256)$ with $\bar \phi = 0$ and 0.3 for $\alpha = 30$, 3000, and 300000. As $\alpha$, which is related to the total chain length of the copolymer, is decreased, the dynamics of the CAC model play a dominant role. These results are consistent with those in ^[18].

3.4.   3D simulation

To observe phase separation in 3D, we use an initial condition as follows:

where $rand(x, y, z)$ is a random number between $-0.001$ and $0.001$ at the grid points. We set $M = 1$, $\epsilon = 0.02$, $\alpha = 300000$, and $\Delta t = \frac{1}{8}$, and evolve $\phi(x, y, z, t)$ for $0 < t \leq 384$. Figure 7 (a) shows the evolution of an isosurface of $\phi(x, y, z, t) = 0$ for $\bar \phi = 0$, 0.26, and 0.4. The final solution forms a gyroidal shape for $\bar \phi = 0$, mixed gyroidal and spherical shapes for $\bar \phi = 0.26$, and a pure spherical shape for $\bar \phi = 0.4$. Figure 7 (b) shows the evolution of $\tilde {\mathcal E}$ for $\bar \phi = 0$, 0.26, and 0.4, which demonstrates that the proposed scheme is energy stable.

4.   Conclusions

We presented the linear, second-order, mass conservative, and energy stable scheme based on the Crank–Nicolson formula. In the scheme, $f(\phi)$ and $\frac{1}{|\Omega|} \int_\Omega f(\phi) \: d{\mathbf{x}}$ were explicitly treated, which made the scheme linear, and the energy stability was guaranteed by adopting the truncated double-well potential that satisfies (1.4) and by adding $A \Delta t \delta_t \phi^{n+1}$ and $B (\delta_t \phi^{n+1} - \delta_t \phi^n)$. We proved that the scheme is mass conservative and energy stable under $A \geq \frac{ML^2}{16}$ and $B \geq \frac{L}{2}$. The numerical results were consistent with the experimental results in ^[17] and confirmed the superior performance of the proposed scheme for phase separation of diblock copolymers.

Appendix

The MATLAB code for the numerical scheme in 2D is provided below.

Use of Generative-AI tools declaration

The author declares they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

The corresponding author (H.G. Lee) thanks the reviewers for the constructive and helpful comments on the revision of this article and was supported by the National Research Foundation of Korea (NRF) grant funded by the Korea government (MSIT) (No. RS-2022-NR069708).

Conflict of interest

The author declares no conflicts of interest in this paper.