An adaptive finite element method based on Superconvergent Cluster Recovery for the Cahn-Hilliard equation

1.   Introduction

The Cahn-Hilliard equation is solved in this study using the adaptive finite element approach.

In the equation above, we use to $\Omega$ signify a limited domain in $R^d(d = 2\ {\rm or} \ 3)$ with Lipschitz boundary $\partial\Omega$, and ${\bf n}$ stands for the unit outward normal to ${\bf \Omega}$, ${\bf u}$ also stands for the phase-field variable that will be solved. $f(u)$ is a nonlinear function. $\varepsilon$ is a constant. $T$ is the ultimate time.

The Cahn-Hilliard (CH) equation, often known as the common phase field model, was first published in 1958 in a key paper by Cahn and Hillard's ^[2], which examined the thermodynamic phenomenon of mutual diffusion of two substances (such as alloys and polymers). In general, this equation can be used to express the intricate phase separation and coarsening phenomena in solids, especially in materials science and fluid dynamics, for examples, see ^[3,4,5,6] and the references therein. The Cahn-Hilliard equation, which holds multiple time scales and spatial scales, is a rigid and nonlinear fourth-order PDE. So it is challenging to identify the precise solution.

The majority of studies in recent years have focused on finite difference techniques or the Fourier-spectral approach (periodic boundary), and many authors have devoted their time to examining different Cahn-Hilliard equation versions, including the viscous Cahn-Hilliard equation, the surface Cahn-Hilliard equation, the CHNS equation, etc. and research issues from space and time that are related. Barrett and Elliott presented to study the Cahn-Hilliard system by using conforming and nonconforming finite element methods ^[7,8,9]. Shen and Xu described in ^[10] the treatment of the Cahn-Hilliard equation by the scalar auxiliary variable (SAV) method. Zhao and Xiao combined the surface finite element method (SFEM) with stabilized semi-implicit model to solve the surface Cahn-Hilliard model ^[11]. There are many parallels between the Allen-Cahn equation and the Cahn-Hilliard equation since they are phase field flows of the same energy in different regions, such as energy decreasing. In ^[12], Shen and Yang analyzed the discrete versions of the Allen-Cahn and Cahn-Hilliard equations. Numerous research has been conducted on the Allen-Cahn equation ^[13,14,15], and Chen and Huang applied the Superconvergent Cluster Recovery (SCR) approach for the Allen-Cahn equation ^[16]. The method was also used in ^[17] to solve the CHNS equation, but the Cahn-Hilliard equation has not yet seen this treatment. In ^[17], the Cahn-Hilliard component is treated as a first-order format combing with a time-space adaptive method. In our paper, we analyze the Cahn-Hilliard equation using a second-order time discrete method and convex splitting for the nonlinear factor $f$, combined with the SCR method.

There are two reasons why we use adaptive methods in this study. First, due to the presence of the tiny parameter $\varepsilon$, it produces the phenomenon of the interfacial layer. When meshing on this interfacial layer, different from other parts, we need fine encryption, and the adaptive strategy can effectively solve this problem. Second, it takes a very long time for the numerical simulation of the Cahn-Hilliard model to reach the steady state, therefore using adaptive mesh creation makes sense in terms of time and money savings. At present, numerous writers have examined the Cahn-Hilliard equation's adaptive technique, including ^[11], where the discrete chemical potential is employed as an error estimator to evaluate the fluctuation in numerical energy. Based on the a posteriori error estimations, the standard mesh refinement procedures were used in ^[18] for automatic mesh refining. In ^[17], the SCR based posteriori error estimators were constructed for space discretization of phase field variable and velocity function. The adaptive time-stepping method was applied to change the time step in ^[19,20].

In this work, we first provide an error estimate for the second-order fully discrete technique suggested in ^[1]. The fully discrete scheme is then introduced with an effective time-space adaptive technique called SCR. We need understand that the SCR approach has the same or higher accuracy while being more effective and affordable than other adaptive techniques. The SCR strategy was originally proposed by Huang and Yi ^[21], and they applied it to the Allen-Cahn system ^[16]. Moreover, for the definition of error estimator, the time discretization error estimator is given by the difference of numerical approximation between adjacent time steps, and the spatial discretization error estimator is defined by gradient, that is, the difference between the reconstructed gradient and the numerical gradient.

The remainder of the essay is structured as follows. For the Cahn-Hilliard equation, we offer a second-order fully discrete scheme in Section 2, after which we examine the error estimation and unconditional energy stability. In Section 3, we introduce the characteristics of the superconvergent cluster recovery operator, and define the error estimation operator in time and space. Several numerical tests illustrate the effectiveness of our adaptive procedure in Section 4.

2.   This energy-stabilized discrete model

First of all, we introduce a new variable ${\bf \mu} = :-\varepsilon\Delta{\bf u}+\frac{1}{\varepsilon}f(u)$ and rewrite Eq (1.1) as follows.

where $\varepsilon > 0$ is a small parameter, which represents the thickness of the transition interface between materials, and ${\bf u}$ means the phase field variable, ${\bf u_{0}}$ $\in$ $H^{1}(\Omega)$ as the initial value at $t = 0$. The solution is driven into the two pure states $u = \pm1$ by the nonlinear part $f(u) = F'(u)$, and $F(u)$ is a double-well potential function, which is given by $F(u) = \frac{1}{4}(u^{2}-1)^{2}$.

Actually, the Cahn-Hilliard equation could be viewed as the $H^{-1}-$gradient flow of the free energy $E(u)$:

Furthermore, by a simple mathematical derivation, for the Cahn-Hilliard equation, the following energy law can be derived:

that the energy-decreasing property is satisfied by the free energy $E(u)$, over time,

Besides, the Cahn-Hilliard equation satisfies the law of conservation of mass ^[2].

Next, we will present a few notations that will be used all across the paper before we get started. The standard Sobolev space will be indicated by the notation $H^m = W^{m, 2}(\Omega)$, and its associated norm of $H^m$ can be represented by $\|\cdot\|_m$. As well as $\|\cdot\|$ and $(\cdot, \cdot)$ respectively stand in for norm and the inner product of $L^2$. Besides, we use the symbol $\|\cdot\|_{\infty}$ to signify the $L^{\infty}$ norm.

Let size $h$ be the largest element diameter for every $e_h$ and $\mathcal {T}_{h} = \{e_{h}\}$ be the set of regularly shaped triangles. Let $V_{h}$ be the corresponding space in finite dimensions for piecewise linear continuous functions:

and its basis functions are the standard Lagrange basis functions $\phi_z (z\in\mathcal {N}_h,$ where the $\mathcal {N}_h$ represents the set of the triangular vertices). Naturally, we denote $V^0_h = V_{h}\bigcap H^1_0(\Omega)$.

In addition, for the real number $p\geq0$ and satisfying the condition $v(t)$, we define the following norm:

Then, by using the above symbol, we first introduce the weak form of Eq (2.1), which is as follow: find ${\bf u}\in L^{\infty}(0, T;H^{1}(\Omega))$, $w\in L^2(0, T;H^1(\Omega))$

To solve Eq (2.6), the semi-discrete system needs to find ${\bf u}_{h} \in C^{1}(0, T;V_{h})$

where ${\bf u}^0_h$ is an approximation of ${\bf u}_0$ in $V_h$.

Divide $[0, T]$ into $N_{st}$ subintervals $I_n: = (t^{n-1}, t^n]$, uniformly, $0 = t^0 < t^1 < ... < t^N = T$. Let $\triangle t$ represent the time step size so $T = N\triangle t$, $n = 1, ..., N$. The fully discrete scheme for the nonlinear approximation, which is developed from a modified Crank-Nicolson approach, is as follows: find $u^n_h \in V^n_h, n = 1, 2, ..., N$

where $\widehat{u^n_h}$, and $f(u^n_h, u^{n-1}_h)$ are defined as

as well as $\mu^{n-\frac{1}{2}}$ is the approximation at the midpoint $t^{n-\frac{1}{2}} = (t^{n-1}+t^n)/2$ (directly computed), and $\Pi_n$ represents the interpolation into the finite element space $V_h$. Next, we consider the mass conservative, taking $w_h = 1$ in (2.8), the following results can be easily obtained

Obviously, the fully discrete scheme is mass conservative.

In addition, We use the convex partitioning method to divide $f$ into two parts, i.e., the linear convex part and the nonlinear concave part. The linear convex part is then handled by the modified Crank-Nicolson method, while the linear convex part is handled by the Midpoint approximation (MP) method. So far, the modified Crank-Nicolson method has been applied several times to the Cahn-Hilliard equation to deal with the nonlinear concave terms, such as ^{[22,23,24,25]}, and the linear convex terms in these works of literature are treated by BDF2.

Lemma 2.1. (Discrete Gronwall lemma) ^[26]. Let $C$ and $\triangle t$ be non-negative constants, and $a_k, b_k, c_k, d_k$ be non-negative sequences satisfying

then

Theorem 2.1. The solution of the fully discrete scheme Eqs (2.8) and (2.9), which is unconditionally energy stable, satisfies

Proof. Taking $w_h = \triangle t\mu^{n-\frac{1}{2}}_h$ and $v_h = u^{n}_h-u^{n-1}_h$ in Eq (2.8) we obtain

where

by integrating the upper formula, we obtain

Also for energy we have

the equation is finally obtained

which means that the energy-decreasing feature is maintained by the provided fully discrete scheme Eqs (2.8) and (2.9). We omit $\Pi_n$ since it implies interpolation into the finite element space $V^n_h$ and has no bearing on the theoretical derivation procedure.

To derive the error estimates, we define the elliptic projection operator as $P_h:H^1\rightarrow V_h$, which satisfies

From the literature ^[27], we know the following two inequalities hold,

Then, we denote the error functions as,

We may quickly derive by using the Taylor expansion with integral residuals and Young's inequality ^[12].

where

Theorem 2.2. Let $u^n_h$ and $u(t^n)$ represent, respectively, the solutions of Eqs (2.8), (2.9) and (2.1). After that, for $u \in C(0, T;H^2(\Omega))$, $u_t\in L^2(0, T;H^2(\Omega))\bigcap L^2(0, T;L^4(\Omega))$, $u_{tt}\in L^2(0, T;H^1(\Omega))$, and $u_{ttt}\in L^2(0, T;H^{-1}(\Omega))$ we get

where

Proof. Subtracting Eq (2.8) from the weak form of Eq (2.6) on $\Omega$ at $t^{n-\frac{1}{2}}$, we obtain

then taking $w_h = \frac{1}{2}\varepsilon\triangle t(\tilde{e}^n+\tilde{e}^{n-1})$ and $v_h = \triangle t\bar{e}^{n-\frac{1}{2}}$ respectively, and summing up the two identities above, we gain

where

the Cauchy inequality with $\varepsilon$ and Young's inequality were used to estimate ${\bf I }$, ${\bf II}, \; {\bf III}$, and ${\bf V}$ as following:

For simplicity, we denote $u(t^n) = u_n$ and $u(t^{n-\frac{1}{2}}) = u_{n-\frac{1}{2}}$. We consider the fourth term Ⅳ as follows:

where

We continue our analysis of these four terms using the Taylor expansion with integral remainders, Cauchy-Schwarz inequality, and Eqs (2.12) and (2.13), we arrive at:

where $\xi_1$ lies between $(u_n+u_{n-1})/2$ and $u_{n-\frac{1}{2}}$, $\xi_2$ lies between $u_n$ and $u_{n-1}$, $\zeta_1$ and $\zeta_2$ are located between $u_{n-1}$ and $u_n$.

where $\xi_3$ lies between $u_n$ and $u_{n-1}$, $\xi_4$ lies between $u_{n-1}$ and $u^{n-1}_h$, and $\xi_5$ lies between $u_n+u_{n-1}$ and $u^n_h+u^{n-1}_h$. To guarantee the existence of $\xi_3, \xi_4, \xi_5$, we need the condition ^[28] (Suppose $C$, a positive constant that is unaffected by $\triangle t$ or $h$, exists).

Further, combining the above inequalities ${\bf IV}_i, i = 1, 2, 3, 4$ into ${\bf IV}$, we arrive at

Moreover, when we combine the terms ${\bf I, II, III, IV}$ and ${\bf V}$, we arrive at the following estimates,

Using Eq (2.11), we may add up the aforementioned inequality for $n = 1, ..., k\; (k\leq T/\triangle t)$, and we find,

Lastly, by using the triangular inequality and the discrete Gronwall lemma to the above inequality, we may conclude the proof.

3.   A posteriori error estimation based SCR

The SCR technique takes derivatives to determine the recovered gradient at recovered locations after fitting a linear polynomial to solution values in a set of suitable sampling points around the vertex. In comparison to other methods, the SCR recovery technique implements adaptive algorithms more simply while saving on computing costs. The crucial component of the SCR technique is the introduction of a posterior error estimate operator. The process for getting the recovered gradient, which is broken down into three parts, is described in the sections that follow.

${\bf Step 1:}$ We plan to restore the gradient of an interior vertex $a = a_0\in \mathcal {N}_h$ in $\Omega$ ($\mathcal {N}_h$ represents the mesh nodes), try to choose some points symmetrically distributed around $z$, and then form a new point set $A = \{a_i = (x_i, y_i)\}$, $1\leq i\leq n$ ($n\geq4$). Specially, if select mesh nodes, we will get better results.

${\bf Step 2:}$ The recovery gradient operator $G_h:V_h\rightarrow V_h\times V_h$ is defined as follows

Let $K_z$ represent the sample points of convex polygons, the linear polynomial $p_a(x, y)$ can be found

In order to overcome the instability caused by small parameter h, we denote $F$ by

where $h: = max\{|x_i-x_0|, \; |y_i-y_0|\}$, $i = 1, 2, ..., n$. Then the fitting polynomial can be written as

with

This $\hat{m}$ is the coefficient vector, it satisfies the linear systems

where

Finally, the recovered gradient can be derived

${\bf Step 3:}$ the recovered gradient $G_hu$ on $\Omega$ is obtained by the interpolation

Next, the following are the properties of the operator $G_h$ was proved in ^[21].

a. For $e_h\in \mathcal {T}_h$, there is a constant $C$ which is independent of the value of $h$, satisfied

where $K = \bigcup^3_{i = 1}\mathcal {K}_i$ with $\mathcal {K}_i$ indicates the components containing the sample points coming from the $i$th vertex of $e_h$.

b. The recovering point $a = a_0 = (x_0, y_0)$ is the center of the circle, and the sampling points $a_i = (x_i, y_i), i = 1, 2, ..., n (n\geq 4)$ equally distributed around it. we derive

We use $\eta^n_t$ and $\eta_{h, e_h}$ as discretization error indicators in time and space respectively, and then adjust the time step and mesh adaption in the algorithm.

During SCR discretization, the time step changes with the error estimate, and the grid size changes according to the time step, so each $t_n$ constructs a new grid, and continues to construct a new space $V_h$ on that grid and $\Pi_n$ represents the interpolation into the finite element space $V_h$.

Remark 3.1. If the gradient is on the boundary $\partial\Omega$, we can treat it in some methods such as a higher order extension, extrapolation, and take the average value.

Remark 3.2. The sample points we selected are placed with maximum symmetry around $a$. On the basis of this, we talk about the sample points in the next two situations. One is that the approximation order can be raised if they are also nodes. The other is that it will only improve recovery accuracy if they are not nodes.

Remark 3.3. The SCR strategy can be used not only for the Crank-Nicolson approximation of time-discrete the Cahn-Hilliard equations, but also for other time-discrete formats of the Cahn-Hilliard equations, for example, we can apply to the BDF approximation mentioned in the literature ^[29][30], and its stability and unique solvability have been demonstrated. Unlike the format in this paper, the format has three levels of time, so we need $u^n_h$ and $u^{n-1}_h$ to approximate $u^{n+1}_h$.

4.   Numerical experiments

To demonstrate the viability of our strategy, some numerical findings are offered in this section. We discuss a few characteristics of the Cahn-Hilliard Eqs (2.8) and (2.9), which are fully discrete and have various beginning value requirements.

Example 4.1. In the first test, we use the initial condition listed below to investigate Eq (2.1).

where the parameter $\Omega = [-1, 1]\times[-1, 1], \, and\; \varepsilon = 0.01$.

Sequences of the mesh adaption results and the numerical solutions for six distinct time steps are shown in Figure 1. We can observe that when the node count decreases, the mesh refinement moves in accordance with the level-set zeros. Figure 2 shows the energy level dropping over time. Furthermore, it is also possible to perceive intuitively how the number of nodes and time steps change over time. We see that during the stage of energy's rapid evolution, a small time-step is chosen to capture the change in the numerical solution, and as the system settles into a stable state, the time steps increase.

Example 4.2. In terms of the Cahn-Hilliard Eq (2.1). Allowing for the following initial condition, let $\Omega = [-1, 1]\times[-1, 1]$.

with the parameter $\varepsilon = 0.01$.

In this illustration, Figure 3 shows a succession of mesh adaptation results and their numerical solutions. It is therefore evident via six different time steps that the mesh transformation adapts the rule of the zeros level-set and the number of nodes is decreasing. We infer from Figure 4 that the energy that decays over time, the number of nodes, and the number of time steps all change with time. This numerical result therefore shows that the time steps will grow greater than before when the system stabilizes.

Example 4.3. In the third example, we apply the following initial value to the Cahn-Hilliard Eq (2.1)

where the parameter $\Omega = [-1, 1]\times[-1, 1], \, and\; \varepsilon = 0.01$.

Figure 5 shows the various adaptive meshes and the associated numerical solutions. It is evident that the mesh transformation applies the level-set rule for zeros. The discrete-time history in Figure 6 also shows that the energy degrades over time. Additionally, Figure 6 shows the connection between time and the quantity of nodes (time-step). In the early stages of the simulation, small steps are used to record changes in the numerical solution, and when the system stabilizes in the latter stages, bigger time steps are used.

Example 4.4. In the last example, Eq (2.1) with the following initial value was analyzed

where the parameter $\Omega = [-1, 1]\times[-1, 1], \, and\; \varepsilon = 0.01$.

In Figure 7, we plot the different adaptive meshes and their contour plots of the numerical solutions. Figure 8, it displays energy is getting more and more small over time and the distribution of the time steps as well as nodes. We come to the conclusion quite similar to the previous examples.

Remark 4.1. When $\varepsilon$ takes a value less than $0.01$, the change of node number throughout the process still follows the law of zero level set, and the energy is still declining. But the time step will become smaller and the computational cost of the program will become larger, so finally, in the paper we choose $\varepsilon = 0.01$ as the experimental result.

Table 1 displays the whole discrete format's spatial and temporal precision Eqs (2.8) and (2.9). We take the initial value $u_{0}(x, y) = 0.01rand(x, y)$, the parameter $\varepsilon = 0.01$, and $T = 1$. Since the precise solution to the initial equation is unknown, the numerical error $e_u^h: = u^h({\bf x}, T)-u^\frac{h}{2}({\bf x}, T)$ will be determined by the difference between the coarse and fine grids, and the convergence order $log_2(\|e_u^h\|/\|e_u^\frac{h}{2}\|)$ will be determined by the ratio of the errors. The table shows that the accuracy converges to 2, which is consistent with our theoretical findings for both spatial and temporal accuracy.

5.   Conclusions

The fundamental idea of this paper is to obtain the spatial discretization operator by SCR method, and then use it as the index of spatial discretization of the Cahn-Hilliard equation. Additionally, the index of temporal discretization is determined by the variation in approximate solutions between adjacent time steps. First of all, in order to solve the Cahn-Hilliard equation, we first develop a second-order scheme that is unconditionally energy stable. Then, the fully discrete system has second-order accuracy both time and space, as determined by error estimation. Finally, these numerical findings are used to demonstrate the efficacy of this approach.

Acknowledgments

Tian's research was supported by Shanxi Scholarship Council of China (No. 2021029) and the 2021 Shanxi Science and Technology Cooperation and Exchange Special Program (No. 202104041101019). Chen's reserach was partially supported by NSFC Project (12201010).

Conflict of interest

No potential conflict of interest was reported by the authors.