An Explicit Adaptive Finite Difference Method for the Cahn-Hilliard Equation.

In this study, we propose an explicit adaptive finite difference method (FDM) for the Cahn-Hilliard (CH) equation which describes the process of phase separation. The CH equation has been successfully utilized to model and simulate diverse field applications such as complex interfacial fluid flows and materials science. To numerically solve the CH equation fast and efficiently, we use the FDM and time-adaptive narrow-band domain. For the adaptive grid, we define a narrow-band domain including the interfacial transition layer of the phase field based on an undivided finite difference and solve the numerical scheme on the narrow-band domain. The proposed numerical scheme is based on an alternating direction explicit (ADE) method. To make the scheme conservative, we apply a mass correction algorithm after each temporal iteration step. To demonstrate the superior performance of the proposed adaptive FDM for the CH equation, we present two- and three-dimensional numerical experiments and compare them with those of other previous methods.

Introduction

The Cahn–Hilliard (CH) equation represents the dynamics of spinodal decomposition that occurs when a single thermodynamic phase spontaneously separates into two phases (Cahn and Hilliard 1958; Grant 1993):where is a phase field at space and time t, is a positive constant, and . The boundary condition is the homogeneous Neumann boundary condition: and for , where is the outer unit normal vector. Because the CH equation can efficiently handle topological changes, it has been widely studied and applied in various applications. However, because the CH equation is a fourth-order nonlinear partial differential equation (PDE), the explicit time discretization requires a very small time step for stability. In order to overcome the strict time step constraints, various semi-implicit (Zhu et al. 1999; Li and Qiao 2017; Li et al. 2016, 2021) and implicit (Beneŝová et al. 2014; Bosch et al. 2014; Xu et al. 2019) schemes have been proposed to solve the CH equation with good numerical stability. Splitting-type numerical schemes have been used for variable studies (Chen et al. 2012; Li et al. 2018; Cheng et al. 2019; Chen et al. 2020; Meng et al. 2020; Hao 2021) with energy stability. Several different methods have been utilized to numerically solve the CH equation. Zhai et al. (2021) proposed a novel high-order linearly operator splitting method for the nonlocal viscous CH equation with the spectral deferred correction method. Their method can avoid the nonlinear iteration and obtain high temporal accuracy. Feng et al. (2020) considered the CH equation with an imposed advection term. They studied the effects of the imposed advection on phase separation of the CH equation. Ainsworth and Mao (2017) investigated the fractional CH equation with fractional free energy. They presented well-posedness and several numerical examples of the fractional CH equation. Li et al. (2016a) used the phase-field model to the minimized CH dynamics. Mohammadi and Dehghan (2019) presented an adaptive time algorithm for the CH equation. Fu and Han (2021) used a finite element method (FEM) to overcome the challenges of solving the CH equation with a degenerate mobility function. Theljani et al. (2020) proposed the image inpainting method which solves the CH equation on an adaptive mesh.Figure 1 shows an adaptive mesh (Theljani et al. 2020).

Fig. 1

a Uniform mesh and b adaptive mesh. Reprinted from (Theljani et al. 2020) with permission from SAGE Publishing

Chen et al. (2016) solved the CH-type diffusion equation and the Navier–Stokes equation using a space-time adaptive finite difference method. Isotopic and strongly anisotropic CH systems were solved using the multigrid method on an adaptive grid (Chen et al. 2018). Figure 2a and b shows the adaptive meshes (Chen et al. 2016, 2018) for two- and three-dimensional spaces, respectively.

Fig. 2

a and b are 2D and 3D adaptive meshes, reprinted from (Chen et al. 2016) and (Chen et al. 2018) with permission from Elsevier, respectively

Compared to the unstructured methods, the adaptive mesh refinement (AMR) framework can provide spatial multi-resolution (Berger and Oliger 1984), which imposes a good choice for the large-scale simulation with challenging physical issues (Li et al. 2016b; Zhou and Xie 2021). There are many applications for the AMR algorithm, i.e., crystal growth (Li and Kim 2012, 2017), thin film (Li et al. 2014; Sun et al. 2007), nonlinear wave (Dohnal and Uecker 2016), and computational fluid dynamics (Berger and Colella 1989). Koliesnikova et al. (2021) proposed a hierarchical mesh refinement scheme driven by a posteriori error estimator. Li et al. (2021) proposed an implicitly discretized surface method with an AMR method. To overcome the difficulties of imposing the boundary conditions by considering the skin depth, Jung and Yoo (2021) used an AMR framework for microwave applications of metallic structures. Grave and Coutinho (2021) used a PDE model to investigate the dynamics of COVID-19 using an AMR method. Their model can represent various spatial scales.

Kay and Welford (2006) presented a multigrid FEM for the CH equation. Banas and Nürnberg (2008) presented an AMR method for the CH equation. Ceniceros and Rom (2007) presented a nonstiff, fully AMR method for the CH equation. Wise et al. (2007) presented a second-order accurate and adaptive FDM to solve the CH equation in 2D and 3D spaces. Stogner et al. (2008) proposed a variational formulation and FEM with an AMR method. The fourth-order spatial accurate compact schemes were developed for the CH equation in 2D (Lee et al. 2014) and 3D (Li et al. 2016c) spaces. The above-mentioned methods are based on a semi-implicit splitting or implicit scheme, in which a multigrid method or other numerical solvers should be performed on coarse and fine grid levels. Therefore, the explicit AMR-based method is much simpler to implement.

The main purpose of this study is to present an explicit adaptive FDM for the CH equation. The primary advantage of the proposed method is its simplicity compared to existing adaptive numerical schemes for the CH equation which use complex data structure and implicit solvers.

The outline of this article is as follows. In Sect. 2, the proposed numerical method is described. In Sect. 3, we present several computational experiments to confirm the superior performance of the proposed explicit AMR method. We conclude this paper in Sect. 4.

Numerical Method

We describe the proposed numerical solution algorithms on the adaptive domain for the 2D and 3D CH equations using the Saul’yev scheme (Yang et al. 2022) and adaptive methodology (Jeong et al. 2021).

Two-Dimensional Algorithm

For the completeness of exposition, we briefly describe the explicit finite difference scheme for the CH equation in the 2D domain . For more details, see (Yang et al. 2022). Let , where the spatial space step is . Here, and are the numbers of the grid points. Let , where is the time step for a nonnegative integer n. First, we discretize Eq. (1) using the linear convex splitting scheme (Li et al. 2017):where and . We define a temporary narrow domain aswhere is a parameter andA space-time adaptive narrow-band domain is defined using :for some positive integer m, which makes buffer points. We can control the narrow-band domain by adjusting and m values appropriately. We define doubly layered exterior boundary points of asFor better understanding of the narrow-band domain construction, let us consider the following circular shape at time :which is shown in Fig. 3a and the red line is the contour of at zero level. Figure 3b, c, and d shows with -plane, temporary domain , and space-time adaptive narrow-band domain (closed circle) with the exterior boundary points (open circle), respectively. We use the values defined at for the exterior boundary values when we solve Eq. (6) on the interior boundary points. Then, using the Saul’yev method (Yang et al. 2022), we have the following updating scheme: For . If , thenwhere . The other 7 ‘for loop’ cases are similarly defined. Next, we adopt the adaptive methodology, which was developed for the Allen–Cahn equation, see (Jeong et al. 2021).

Fig. 3

Schematic illustration for , and : a mesh plot of with the zero level contour (red line), b and -plane, c temporary domain , and d space-time adaptive narrow-band domain (closed circle) with the exterior boundary points (open circle) (Color figure online)

Because the Saul’yev-type scheme is generally not conservative, we adopt a mass correction step (Jeong et al. 2020):We should note that the mass correction is only done on the narrow-band domain because the values of outside the narrow-band domain are unchanged.

Remark 1

Solving problems using the FDM in adaptive meshes is very complex and difficult; nevertheless, there are studies that have presented new numerical ideas and convincing numerical results on this interesting scientific problem (Feng et al. 2018; Wise et al. 2011). However, because most adaptive methods have very complex relational expressions, numerical methods have been verified through numerical experiments rather than analytically verified.

Three-Dimensional Algorithm

The 3D numerical solution algorithm and adaptive grid generation are straightforward extension of the 2D ones. Let . Let be the discrete domain, where . Here, , and are the numbers of the grid points. Let . We start with the following method (Li et al. 2017):where and . Then, using the Saul’yev method (Yang et al. 2022), we have the following updating scheme:where . The other 47 ‘for loop’ cases are similarly defined. A space-time adaptive narrow-band domain is defined as whereFigure 4 illustrates a schematic of narrow-band domain on 3D space. It shows the temporal evolution of (dots) with isosurface of solutions at zero level.

Fig. 4

a–c Snapshots of the temporal evolution of narrow-band domain (dots) in 3D space

To make the numerical scheme conservative, we apply the mass correction algorithm after each temporal iteration:In Eqs. (2) and (8), we only considered a first-order linear convex splitting scheme. It can be extended to second-order linear convex splitting schemes (Guo et al. 2016, 2021; Chen et al. 2019; Dong et al. 2020). For example, we may use the following second-order method (Guo et al. 2016): for given ,where

Computational Tests

We perform various computational simulations using the proposed numerical solution algorithm on 2D and 3D spaces. Through the computational experiments, we demonstrate the efficiency and accuracy of our proposed numerical solution algorithm. For numerical simulations, we define the interface layer parameter as follows (Jeong et al. 2021):where l is a positive integer. Unless otherwise stated, numerical experiments are performed in () with grid points , uniform spatial step , and .

Stability Test

In this section, we test the time step stability of the proposed adaptive method. In 2D and 3D spaces, we consider square and cubic shapes as initial conditions, respectively. We define a maximum time step , which is a numerical maximum time step that guarantees the non-blow-up of the numerical solutions. The condition under which the numerical solution blows up is defined when , where I is the index set of grid points. In 2D and 3D computational domains , where d is dimension, the given initial conditions are as follows:andrespectively. Tables 1 and 2 list the maximum time step in 2D and 3D space, respectively, with three different spatial grid sizes and 1/256. The stability test confirms that the proposed adaptive method can use sufficiently large time steps for the stable numerical solutions.

Table 1

Stability test results in 2D space

Case	Δt_max
h=1/64	8.8199e-4 ≈14790h^4
h=1/128	2.2140e-4 ≈59432h^4
h=1/256	5.6189e-5 ≈241330h^4

Table 2

Stability test results in 3D space

Case	Δt_max
h=1/64	1.1704e-4 ≈1964h^4
h=1/128	6.3239e-5 ≈16976h^4
h=1/256	1.5285e-5 ≈65650h^4

Two-Dimensional Space

In 2D space, we compare solutions of the entire domain with those of adaptive domains and investigate effects of the parameters m, , and K which generate the adaptive domain.

In Fig. 5, we observe the temporal evolution with the following square shape initial condition,

Fig. 5

Temporal evolution of the contours of at zero level with square initial condition. a–c are snapshots of solution for entire and adaptive domains at , , and , respectively. Here, , , and are used

We compare the solutions of the CH equation computed on the entire domain (solid line) with the adaptive domain (solid line with markers). The solutions on the adaptive domain fit well with the solutions on the entire domains. For efficient computations, buffer size m which adjusts the narrow-band domain should be appropriately controlled. To find a proper the buffer size m, we investigate the evolution of the solutions on the adaptive narrow-band domain according to the value of m. In Fig. 6, we experiment with a square spiral shape. As a result of Fig. 6, is sufficient to accurately represent the solution of the CH equation.

Fig. 6

Snapshots of solution and narrow-band domain for a–b adaptive domain and , c entire domain. d Overlapped contour of solutions at zero level. Left to right column, , , and . Here, and are used

To compare CPU times for the entire domains with those for adaptive domains, we define a period K of updating the adaptive domain at every K time steps. That is, we update the computational domain when for , where K is a positive integer. In Fig. 7, we show the comparison results for entire and adaptive domains. The CH equation is solved on a rectangular domain with parameter values: , , , , and . We investigate the dynamics until the solution of the CH equation becomes circular from the rectangular initial condition. Figure 7a shows CPU times for entire and adaptive domains. It shows that as K increases, the CPU times decrease. Figure 7b shows contours of at zero for the entire and adaptive domains with the rectangular initial condition. Solutions in the adaptive domains are similar to the solution in the entire domain up to . As shown in Figs. 7a and b, it seems reasonable to adopt for efficient and accurate computation.

Fig. 7

On a computational domain , for a rectangular initial condition, a CPU times and b contours of at zero for entire domain and adaptive domains with , and 10000. Here, , , , and are used

Three-Dimensional Space

We perform several computational experiments in 3D space. In 3D space, for stability and efficient computation, we use , , , and . Figure 8 displays the temporal evolution of the cubic shape initial condition as follows:

Fig. 8

Temporal evolution isosurface of at zero level with cube initial condition in 3D space. a–c are snapshots of solution for adaptive domain at , , and , respectively

To examine the efficiency of adaptive narrow-domain, a shape with irregularity is considered as the initial condition as shown in Fig. 9. In the case of the shape of Fig. 9, higher efficiencies can be obtained when using the adaptive narrow-band domain, compared to the case of solving the problems on the entire domain. Let us define the 3D discrete energy functional in the discretized space domain aswhere

Fig. 9

Temporal evolution of the discrete free energy functional and dynamics with respect to the number of time iterations. Here, the solid line represents the energy curve, open circles represent the specific moments when , , and . From left to right, the figures inserted into the graph are snapshots of the zero level isosurface of solution with the irregular initial condition

In Fig. 9, the graph shows the temporal evolution of discrete free energy from to and the figures are snapshots of the zero level isosurface of with the irregular initial condition. We confirm that the discrete free energy functional decreases monotonically as time goes on.

Then, we investigate the robustness of the adaptive narrow-band domain for a shape with complex interface in 3D space. We consider ‘Genus 6 3D surface’ (Ceh Jan) as an initial condition, where the inside of the surface is 1 and the outside is -1, as shown in Fig. 10a. Figure 10a–d shows the temporal evolution of isosurface of at zero level with Genus 6 shape. As shown in Fig. 10, it can be seen that the adaptive narrow-band domain operates well for the CH equation for shapes with complex interfaces such as Genus 6 3D surface.

Fig. 10

Temporal evolution of isosurface of at zero level with Genus 6 shape for initial condition. a–d are snapshot of solution for the adaptive domain at , , , and , respectively

To validate the adaptive procedure for disjoint drops in 3D space, we consider the four separate drops, (see, e.q., Cheng et al. (2018)) consisting of three ellipsoids and one sphere. The initial condition is given asThe isosurface of initial condition at is illustrated in Fig. 11a. From the results in Fig. 11, four disjoint drops merge and eventually evolve into a sphere. The solution obtained by the proposed explicit adaptive method follows the dynamics of the CH equation well even if it is initially irregular or separated.

Fig. 11

a–f Snapshots of the isosurface of numerical solution at zero level with initially four drops, using the adaptive domain with at , , , , , and , respectively

Discussion and Conclusion

In this section, we discuss a relation between the interface layer parameter and the gradient criterion . Then, we draw conclusions of this paper. An equilibrium profile is on the infinite domain, where the phase-field varies from to 0.9 at a distance of about . Therefore, we use Eq. (17) for a transition layer width of approximately lh. We differentiate the equilibrium solution to define the gradient criterion value .Thus, we define the relation between and as follows:which implies that our algorithm is set to recognize the transition layer more accurately. When , using the relation (23), . In our previous numerical tests, we use for simplicity.

In this paper, we presented an explicit adaptive FDM for the CH equation which describes the process of phase separation. To numerically solve the CH equation fast and efficiently, we used the FDM and time-adaptive narrow-band domain. For the adaptive grid, we define a narrow-band domain including the interfacial transition layer of the phase field and solve the numerical scheme on the narrow-band domain which is computed using undivided finite difference. We used the alternating direction explicit method. To make the scheme conservative, we apply a mass correction algorithm after each temporal iteration step. To demonstrate the superior performance of the proposed method, we presented 2D and 3D numerical experiments and compared them with those of the other previous methods. There exist several interesting future extensions of the proposed adaptive approach of the CH equation to other phase field models such as two-phase flows (Guo et al. 2022), multiphase tumor growth (Chen et al. 2014), topology optimization (Bartels et al. 2021; Yu et al. 2021), N-component CH system (Li et al. 2016b), N-component fluid flows (Kim 2009; Xia et al. 2022), curve and surface smoothing (Choi et al. 2017), and dendritic crystal growth (Zhang et al. 2019).