
In this study, we construct an error estimate for a fully discrete finite element scheme that satisfies the criteria of unconditional energy stability, as suggested in [
Citation: Wenyan Tian, Yaoyao Chen, Zhaoxia Meng, Hongen Jia. An adaptive finite element method based on Superconvergent Cluster Recovery for the Cahn-Hilliard equation[J]. Electronic Research Archive, 2023, 31(3): 1323-1343. doi: 10.3934/era.2023068
[1] | Hao Wu . A review on the Cahn–Hilliard equation: classical results and recent advances in dynamic boundary conditions. Electronic Research Archive, 2022, 30(8): 2788-2832. doi: 10.3934/era.2022143 |
[2] | Liupeng Wang, Yunqing Huang . Error estimates for second-order SAV finite element method to phase field crystal model. Electronic Research Archive, 2021, 29(1): 1735-1752. doi: 10.3934/era.2020089 |
[3] | Hao Wang, Wei Yang, Yunqing Huang . An adaptive edge finite element method for the Maxwell's equations in metamaterials. Electronic Research Archive, 2020, 28(2): 961-976. doi: 10.3934/era.2020051 |
[4] | Shan Jiang, Li Liang, Meiling Sun, Fang Su . Uniform high-order convergence of multiscale finite element computation on a graded recursion for singular perturbation. Electronic Research Archive, 2020, 28(2): 935-949. doi: 10.3934/era.2020049 |
[5] | Zengyan Zhang, Yuezheng Gong, Jia Zhao . A remark on the invariant energy quadratization (IEQ) method for preserving the original energy dissipation laws. Electronic Research Archive, 2022, 30(2): 701-714. doi: 10.3934/era.2022037 |
[6] | Suayip Toprakseven, Seza Dinibutun . A weak Galerkin finite element method for parabolic singularly perturbed convection-diffusion equations on layer-adapted meshes. Electronic Research Archive, 2024, 32(8): 5033-5066. doi: 10.3934/era.2024232 |
[7] | Changling Xu, Huilai Li . Two-grid methods of finite element approximation for parabolic integro-differential optimal control problems. Electronic Research Archive, 2023, 31(8): 4818-4842. doi: 10.3934/era.2023247 |
[8] |
Zuliang Lu, Fei Huang, Xiankui Wu, Lin Li, Shang Liu .
Convergence and quasi-optimality of |
[9] | Victor Ginting . An adjoint-based a posteriori analysis of numerical approximation of Richards equation. Electronic Research Archive, 2021, 29(5): 3405-3427. doi: 10.3934/era.2021045 |
[10] | Hsueh-Chen Lee, Hyesuk Lee . An a posteriori error estimator based on least-squares finite element solutions for viscoelastic fluid flows. Electronic Research Archive, 2021, 29(4): 2755-2770. doi: 10.3934/era.2021012 |
In this study, we construct an error estimate for a fully discrete finite element scheme that satisfies the criteria of unconditional energy stability, as suggested in [
The Cahn-Hilliard equation is solved in this study using the adaptive finite element approach.
{ut=−εΔ2u+1εΔf(u)in Ω×(0,T],∂nu=∂n(−εΔu+1εf(u))=0,on ∂Ω×(0,T],u(x,0)=u0in Ω. | (1.1) |
In the equation above, we use to Ω signify a limited domain in Rd(d=2 or 3) with Lipschitz boundary ∂Ω, and n stands for the unit outward normal to Ω, u also stands for the phase-field variable that will be solved. f(u) is a nonlinear function. ε 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 ε, 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.
First of all, we introduce a new variable μ=:−εΔu+1εf(u) and rewrite Eq (1.1) as follows.
{ut=Δμin Ω×(0,T],μ=−εΔu+1εf(u)in Ω×(0,T],∂nu=∂nμ=0on ∂Ω×(0,T],u(x,0)=u0in Ω, | (2.1) |
where ε>0 is a small parameter, which represents the thickness of the transition interface between materials, and u means the phase field variable, u0 ∈ H1(Ω) as the initial value at t=0. The solution is driven into the two pure states u=±1 by the nonlinear part f(u)=F′(u), and F(u) is a double-well potential function, which is given by F(u)=14(u2−1)2.
Actually, the Cahn-Hilliard equation could be viewed as the H−1−gradient flow of the free energy E(u):
E(u)=∫Ω(1εF(u)+ε2|∇u|2)dx, | (2.2) |
Furthermore, by a simple mathematical derivation, for the Cahn-Hilliard equation, the following energy law can be derived:
ddtE(u)=−∫Ω|∇(−εΔu+1εf(u))|2dx≤0, | (2.3) |
that the energy-decreasing property is satisfied by the free energy E(u), over time,
E(u)(tn+1)≤E(u)(tn),∀n∈N. | (2.4) |
Besides, the Cahn-Hilliard equation satisfies the law of conservation of mass [2].
ddt∫Ωu(x,t)dx=0,0≤t≤T. | (2.5) |
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 Hm=Wm,2(Ω), and its associated norm of Hm can be represented by ‖⋅‖m. As well as ‖⋅‖ and (⋅,⋅) respectively stand in for norm and the inner product of L2. Besides, we use the symbol ‖⋅‖∞ to signify the L∞ norm.
Let size h be the largest element diameter for every eh and Th={eh} be the set of regularly shaped triangles. Let Vh be the corresponding space in finite dimensions for piecewise linear continuous functions:
Vh={v∈H1(Ω):v|eh∈P1(eh),∀eh∈Th}, |
and its basis functions are the standard Lagrange basis functions ϕz(z∈Nh, where the Nh represents the set of the triangular vertices). Naturally, we denote V0h=Vh⋂H10(Ω).
In addition, for the real number p≥0 and satisfying the condition v(t), we define the following norm:
‖v‖Lp(0,T;X)=(∫T0‖v(t)‖pXdt)12. |
Then, by using the above symbol, we first introduce the weak form of Eq (2.1), which is as follow: find u∈L∞(0,T;H1(Ω)), w∈L2(0,T;H1(Ω))
{(ut,w)+(∇μ,∇w)=0∀w∈H1(Ω),(μ,v)=ε(∇u,∇v)+1ε(f(u),v)∀v∈H1(Ω),(u(0)−u0,q)=0∀q∈H1(Ω). | (2.6) |
To solve Eq (2.6), the semi-discrete system needs to find uh∈C1(0,T;Vh)
{(uh,t,wh)+(∇μh,∇wh)=0∀wh∈Vh,(μh,vh)=ε(∇uh,∇vh)+1ε(f(uh),vh)∀vh∈Vh,(uh(0)−u0h0,qh)=0∀qh∈Vh, | (2.7) |
where u0h is an approximation of u0 in Vh.
Divide [0,T] into Nst subintervals In:=(tn−1,tn], uniformly, 0=t0<t1<...<tN=T. Let △t represent the time step size so T=N△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 unh∈Vnh,n=1,2,...,N
{(unh−Πnun−1h△t,wh)+(∇μn−12h,∇wh)=0∀wh∈Vnh,(μn−12h,vh)=ε(∇^unh,∇vh)+1ε(f(unh,un−1h),vh)∀vh∈Vnh,(uh(0)−u0h0,qh)=0∀qh∈V0h, | (2.8) |
where ^unh, and f(unh,un−1h) are defined as
^unh=unh+Πnun−1h2,f(unh,un−1h)=(unh)3+(unh)2Πnun−1h+unh(Πnun−1h)2+(Πnun−1h)34−unh+Πnun−1h2, | (2.9) |
as well as μn−12 is the approximation at the midpoint tn−12=(tn−1+tn)/2 (directly computed), and Πn represents the interpolation into the finite element space Vh. Next, we consider the mass conservative, taking wh=1 in (2.8), the following results can be easily obtained
∫Ωun=∫Ωun−1=∫Ωu0. |
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 △t be non-negative constants, and ak,bk,ck,dk be non-negative sequences satisfying
ak+△tn∑k=0bk≤△tn−1∑k=0dkak+△tn−1∑k=0ck+C0,∀n≥1, |
then
an+△tn∑k=0bk≤exp(△tn−1∑k=0dk)(△tn−1∑k=0ck+C0),∀n≥1. |
Theorem 2.1. The solution of the fully discrete scheme Eqs (2.8) and (2.9), which is unconditionally energy stable, satisfies
E(unh)≤E(un−1h),∀n≥1. |
Proof. Taking wh=△tμn−12h and vh=unh−un−1h in Eq (2.8) we obtain
{(unh−un−1h,μn−12h)+△t‖∇μn−12h‖2=0,(μn−12h,unh−un−1h)=ε2(‖∇unh‖2−‖∇un−1h‖2)+1ε(f(unh,un−1h),unh−un−1h), |
where
1ε(f(unh,un−1h),unh−un−1h)=1ε((unh)3+(unh)2un−1h+unh(un−1h)2+(un−1h)34−unh+un−1h2,unh−un−1h)=1ε∫14(unh)4−14(un−1h)4−12(unh)2+12(un−1h)2dx, |
by integrating the upper formula, we obtain
0=△t‖∇μn−12h‖2+ε2(‖∇unh‖2−‖∇un−1h‖2)+1ε∫14(unh)4−14(un−1h)4−12(unh)2+12(un−1h)2dx. |
Also for energy we have
E(unh)−E(un−1h)=∫1εF(unh)+ε2|∇unh|2−1εF(un−1h)+ε2|∇un−1h|2dx=1ε∫14(unh)4−14(un−1h)4−12(unh)2+12(un−1h)2dx+ε2(‖∇unh‖2−‖un−1h‖2), |
the equation is finally obtained
△t‖∇μn−12h‖2+E(unh)−E(un−1h)=0, |
which means that the energy-decreasing feature is maintained by the provided fully discrete scheme Eqs (2.8) and (2.9). We omit Πn since it implies interpolation into the finite element space Vnh and has no bearing on the theoretical derivation procedure.
To derive the error estimates, we define the elliptic projection operator as Ph:H1→Vh, which satisfies
(∇(Phv−v),∇αh)=0,∀αh∈Vh. | (2.10) |
From the literature [27], we know the following two inequalities hold,
‖v−Phv‖≤Ch2‖v‖2,‖(v−Phv)t‖≤Ch2‖vt‖2. | (2.11) |
Then, we denote the error functions as,
˜en=Phu(tn)−unh,ˆen=u(tn)−Phu(tn),ˉen−12=Phμ(tn−12)−μn−12h,ˇen−12=μ(tn−12)−Phμ(tn−12). |
We may quickly derive by using the Taylor expansion with integral residuals and Young's inequality [12].
‖Rn1‖2s≤△t3∫tntn−1‖uttt(t)‖2sdt,s=−1,0, | (2.12) |
‖Rn2‖2s≤△t3∫tntn−1‖utt(t)‖2s+2dt,s=−1,0, | (2.13) |
where
Rn1=u(tn)−u(tn−1)△t−ut(tn−12), |
Rn2=−△(u(tn)+u(tn−1)2−u(tn−12)). |
Theorem 2.2. Let unh and u(tn) represent, respectively, the solutions of Eqs (2.8), (2.9) and (2.1). After that, for u∈C(0,T;H2(Ω)), ut∈L2(0,T;H2(Ω))⋂L2(0,T;L4(Ω)), utt∈L2(0,T;H1(Ω)), and uttt∈L2(0,T;H−1(Ω)) we get
‖ukh−u(tk)‖+(△tεk∑n=0‖μn−12h−μ(tn−12)‖2)12⩽C(ε,T)(K1(ε,u)△t2+K2(ε,u)h2), |
where
C(ε,T)∼exp(T∖ε),K1(ε,u)=√ε(‖uttt‖L2(0,T;H−1)+‖utt‖L2(0,T;H1))+1√ε(‖utt‖L2(0,T;L2)+‖ut‖L2(0,T;L4)),K2(ε,u)=‖u0‖2+√ε‖ut‖L2(0,T;H2)+1√ε‖μ‖C(0,T;H2)+‖u‖C(0,T;H2). |
Proof. Subtracting Eq (2.8) from the weak form of Eq (2.6) on Ω at tn−12, we obtain
1△t(˜en−˜en−1,wh)+(∇ˉen−12,∇wh)=(Rn1−1△t(I−Ph)(u(tn)−u(tn−1)),wh),(ˉen−12+ˇen−12,vh)=ε2(∇˜en+∇˜en−1,vh)+ε(△Rn2,vh)+1ε(f(u(tn−12))−f(unh,un−1h),vh), |
then taking wh=12ε△t(˜en+˜en−1) and vh=△tˉen−12 respectively, and summing up the two identities above, we gain
‖˜en‖2−‖˜en−1‖2+2△tε‖ˉen−12‖2=△t(Rn1,˜en+˜en−1)−(I−Ph)(u(tn)−u(tn−1),˜en+˜en−1)−2△t(Rn2,ˉen−12)+2△tε2(f(u(tn−12))−f(unh,unn−1),ˉen−12)−2△tε(ˇen−12,ˉen−12):=I+II+III+IV+V, |
where
I:=△t(Rn1,˜en+˜en−1),II:=−(I−Ph)(u(tn)−u(tn−1),˜en+˜en−1),III:=−2△t(Rn2,ˉen−12),IV:=2△tε2(f(u(tn−12))−f(unh,un−1h),ˉen−12),V:=−2△tε(ˇen−12,ˉen−12). |
the Cauchy inequality with ε and Young's inequality were used to estimate I, II,III, and V as following:
I≤△t‖Rn1‖‖˜en+˜en−1‖≤△t2(ε‖Rn1‖2+1ε‖˜en+˜en−1‖2)≤ε△t42∫tntn−1‖uttt(t)‖2−1dt+△t2ε‖˜en‖2+△t2ε‖˜en−1‖2,II≤‖(I−Ph)(u(tn)−u(tn−1))‖‖˜en+˜en−1‖≤ε2∫tntn−1‖(I−Ph)ut(t)‖2dt+12ε‖˜en‖2+12ε‖˜en−1‖2,III≤3△tε‖Rn2‖2+△t3ε‖ˉen−12‖2≤3ε△t4∫tntn−1‖utt(t)‖21dt+△t3ε‖ˉen−12‖2,V≤2△tε(√3‖ˇen−12‖)(1√3‖ˉen−12‖)≤3△tε‖ˇen−12‖2+△t3ε‖ˉen−12‖2. |
For simplicity, we denote u(tn)=un and u(tn−12)=un−12. We consider the fourth term Ⅳ as follows:
IV=2△tε2(f(un−12)−f(unh−un−1h),ˉen−12)=2△tε2(unh−un−1h2−un−12+(un−12)3−(unh)3+(unh)2un−1hh+unh(un−1h)2+(un−1h)34)=2△tε2(unh−un−1h2−un−12+(un−12)3−(un−un−12)3+(un−un−12)3−gn+gn−(unh)3+(unh)2un−1h+unh(un−1h)2+(un−1h)34,ˉen−12):=2△tε2(IV1+IV2+IV3+IV4,ˉen−12), |
where
gn=(un)3+(un)2un−1+un(un−1)2+(un−1)34,IV1:=unh+un−1h2,IV2:=(un−12)3−(un+un−12)3,IV3:=(un+un−12)3−gn,IV4:=gn−(unh)3+(unh)2un−1h+unh(un−1h)2+(un−1h)34. |
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:
‖IV1‖=‖−˜en−ˆen−˜en−1−ˆen−1+un+un−12−un−12‖≤‖˜en+˜en−12‖+‖ˆen+ˆen−12‖+‖un+un−12−un−12‖≤‖˜en+˜en−12‖+‖ˆen+ˆen−12‖+△t3∫tntn−1‖utt(t)‖2dt, |
‖IV2‖=‖(un−12)3−(un+un−12)3‖≤‖3ξ21(un+un−12−un−12)‖≤3ξ21△t3∫tntn−1‖utt(t)‖2dt, |
‖IV3‖=‖(un)3−(un)2un−1−un(un−1)2+(un−1)38‖=‖((un)2−(un−1)2)(un−un−1)8‖≤‖ξ24(un−un−1)2‖≤ξ24△t3∫tntn−1‖u2t(t)‖2dt, |
where ξ1 lies between (un+un−1)/2 and un−12, ξ2 lies between un and un−1, ζ1 and ζ2 are located between un−1 and un.
‖IV4‖=14‖23(un)3+13(un+un−1)3+23(un−1)3−23(unh)3+13(unh+un−1h)3+23(un−1h)3‖=‖16((un)3−(unh)3)+16((un−1)3−(un−1h)3)+112((un+un−1)3−(unh+un−1h)3)‖≤12‖ξ23(un−unh)‖+12‖ξ24(un−1−un−1h)‖+14‖ξ25(un−unh+un−1−un−1h)‖≤C(‖un−1−un−1h‖+‖un−unh‖)≤C(‖ˆen−1+˜en−1‖+‖ˆen+˜en‖)≤C(‖ˆen−1‖+‖˜en−1‖+‖ˆen‖+‖˜en‖), |
where ξ3 lies between un and un−1, ξ4 lies between un−1 and un−1h, and ξ5 lies between un+un−1 and unh+un−1h. To guarantee the existence of ξ3,ξ4,ξ5, we need the condition [28] (Suppose C, a positive constant that is unaffected by △t or h, exists).
‖unh‖∞≤C,∀1≤n≤N. |
Further, combining the above inequalities IVi,i=1,2,3,4 into IV, we arrive at
IV=2△tε2(IV1+IV2+IV3+IV4,ˉen−12)≤2△tε2‖IV1+IV2+IV3+IV4‖‖ˉen−12‖≤2△tε2(‖IV1‖+‖IV2‖+‖IV3‖+‖IV4‖)‖ˉen−12‖≤2△tε2(‖˜en+˜en−12‖+‖ˆen+ˆen−12‖+△t3∫tntn−1‖utt(t)‖2dt+3ξ21△t3∫tntn−1‖utt(t)‖2dt+ξ24△t3∫tntn−1‖u2t(t)‖2dt+C(‖ˆen−1‖+‖˜en−1‖+‖ˆen‖+‖˜en‖))‖ˉen−12‖≤3△tε(‖˜en+˜en−12‖+‖ˆen+ˆen−12‖+△t3∫tntn−1‖utt(t)‖2dt+3ξ21△t3∫tntn−1‖utt(t)‖2dt+ξ24△t3∫tntn−1‖u2t(t)‖2dt+C(‖ˆen−1‖+‖˜en−1‖+‖ˆen‖+‖˜en‖))2+△t3ε‖ˉen−12‖2≤C△t4ε∫tntn−1‖utt(t)‖2dt+C△t4ε∫tntn−1‖u2t(t)‖2dt+C(‖ˆen−1‖2+‖˜en−1‖2+‖ˆen‖2+‖˜en‖2)+△t3ε‖ˉen−12‖2. |
Moreover, when we combine the terms I,II,III,IV and V, we arrive at the following estimates,
‖˜en‖2−‖˜en−1‖2+2△tε‖ˉen−12‖2≤ε△t42∫tntn−1‖uttt(t)‖2−1dt+ε2∫tntn−1‖(I−Ph)utt(t)‖2dt+3△tε‖ˇen−12‖2+△t3ε‖ˉen−12‖2+C△t4ε∫tntn−1‖utt(t)‖2dt+C△t4ε∫tntn−1‖u2t(t)‖2dt+3ε△t4∫tntn−1‖utt(t)‖21dt+△t3ε‖ˉen−12‖2+C(‖ˆen−1‖2+‖˜en−1‖2+‖ˆen‖2+‖˜en‖2)+△t3ε‖ˉen−12‖2. |
Using Eq (2.11), we may add up the aforementioned inequality for n=1,...,k(k≤T/△t), and we find,
‖˜ek‖2−‖˜e0‖2+△tεk∑n=1‖ˉen−12‖2≤12△t4ε‖uttt(t)‖2L2(0,T;H−1)+ε2Ch4‖ut(t)‖2L2(0,T;H2)+3△tεk∑n=1‖ˇen−12‖2+C△t4ε‖utt(t)‖2L2(0,T;L2)+C△t4ε‖ut(t)‖L2(0,T;L4)+3ε△t4‖utt(t)‖2L2(0,T;H1)+Ck∑n=1(‖ˆen−1‖2+‖˜en−1‖2+‖ˆen‖2+‖˜en‖2). |
Lastly, by using the triangular inequality and the discrete Gronwall lemma to the above inequality, we may conclude the proof.
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.
Step1: We plan to restore the gradient of an interior vertex a=a0∈Nh in Ω (Nh represents the mesh nodes), try to choose some points symmetrically distributed around z, and then form a new point set A={ai=(xi,yi)}, 1≤i≤n (n≥4). Specially, if select mesh nodes, we will get better results.
Step2: The recovery gradient operator Gh:Vh→Vh×Vh is defined as follows
(Ghuh)(a)=∇pa(x,y). |
Let Kz represent the sample points of convex polygons, the linear polynomial pa(x,y) can be found
pa(x,y):=argminp1∈P1n∑i=0(p1(ai)−uh(ai))2, |
In order to overcome the instability caused by small parameter h, we denote F by
F:(x,y)→(ψ,φ)=(x,y)−(x0,y0)h, |
where h:=max{|xi−x0|,|yi−y0|}, i=1,2,...,n. Then the fitting polynomial can be written as
pa(x,y)=PTm=ˆPTˆm, |
with
PT=(1,x,y),ˆPT=(2,ψ,φ),mT=(m1,m2,m3),ˆmT=(ˆm1,ˆm2,ˆm3). |
This ˆm is the coefficient vector, it satisfies the linear systems
ATAˆm=ATu, |
where
A=(1ψ0φ01ψ1φ1⋮⋮⋮1ψnφn)andu=(u(a0)u(a1)⋮u(an)). |
Finally, the recovered gradient can be derived
Ghu=∇pa=(m2m3)=1h(ˆm2ˆm3) |
Step3: the recovered gradient Ghu on Ω is obtained by the interpolation
Ghu=∑a∈NhGhu(a)ϕa. |
Next, the following are the properties of the operator Gh was proved in [21].
a. For eh∈Th, there is a constant C which is independent of the value of h, satisfied
‖Ghv‖0,eh≤C|v|1,K∀v∈Vh, |
where K=⋃3i=1Ki with Ki indicates the components containing the sample points coming from the ith vertex of eh.
b. The recovering point a=a0=(x0,y0) is the center of the circle, and the sampling points ai=(xi,yi),i=1,2,...,n(n≥4) equally distributed around it. we derive
|∇v(a)−Ghv(a)|≤Ch2,∀v∈W3∞(eh). |
We use ηnt and ηh,eh as discretization error indicators in time and space respectively, and then adjust the time step and mesh adaption in the algorithm.
ηnt:=‖unh−Πnun−1h‖,ηh,eh:=‖Ghuh−∇uh‖0,eh,η2h=∑eh∈Thη2h,eh. |
During SCR discretization, the time step changes with the error estimate, and the grid size changes according to the time step, so each tn constructs a new grid, and continues to construct a new space Vh on that grid and Πn represents the interpolation into the finite element space Vh.
Remark 3.1. If the gradient is on the boundary ∂Ω, 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 unh and un−1h to approximate un+1h.
{(3un+1h−4Πnunh+Πnun−1h2△t,vh)+(∇wn+1h,∇vh)=0,∀vh∈Vh,(wn+1h,φh)=ε2(∇un+1h,∇φh)+((un+1h)3−2Πnunh+Πnun−1h,φh)+A△t(∇(un+1h−Πnunh),∇φh),∀φh∈Vh, |
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).
u0(x,y)=tanh(((x−0.3)2+y2−0.22)/ε)tanh(((x+0.3)2+y2−0.22)/ε)×tanh((x2+(y−0.3)2−0.22)/ε)tanh((x2+(y+0.3)2−0.22)/ε), | (4.1) |
where the parameter Ω=[−1,1]×[−1,1],andε=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 Ω=[−1,1]×[−1,1].
u0(x,y)=tanh((x2+y2−0.152)/ε)×tanh(((x−0.31)2+y2−0.152)/ε)×tanh(((x+0.31)2+y2−0.152)/ε)×tanh((x2+(y−0.31)2−0.152)/ε)×tanh((x2+(y+0.31)2−0.152)/ε)×tanh(((x−0.31)2+(y−0.31)2−0.152)/ε)×tanh(((x−0.31)2+(y+0.31)2−0.152)/ε)×tanh(((x+0.31)2+(y−0.31)2−0.152)/ε)×tanh(((x+0.31)2+(y+0.31)2−0.152)/ε), | (4.2) |
with the parameter ε=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)
u0(x,y)=tanh(((x−0.3)2+y2−0.252)/ε)tanh(((x+0.3)2+y2−0.32)/ε), | (4.3) |
where the parameter Ω=[−1,1]×[−1,1],andε=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
u0(x,y)=0.01rand(x,y), | (4.4) |
where the parameter Ω=[−1,1]×[−1,1],andε=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 ε 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 ε=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 u0(x,y)=0.01rand(x,y), the parameter ε=0.01, and T=1. Since the precise solution to the initial equation is unknown, the numerical error ehu:=uh(x,T)−uh2(x,T) will be determined by the difference between the coarse and fine grids, and the convergence order log2(‖ehu‖/‖eh2u‖) 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.
h | ‖umh−umh2‖0 | Rate | △t | Rate |
116 | 0.052470400 | 116 | ||
132 | 0.018185700 | 1.528700 | 132 | 2.18632 |
164 | 0.003449260 | 2.398450 | 164 | 2.52708 |
1128 | 0.000945365 | 1.867340 | 1128 | 1.78946 |
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.
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).
No potential conflict of interest was reported by the authors.
[1] |
F. Guillen-Gonzalez, G. Tierra, Second order schemes and time-step adaptivity for Allen-Cahn and Cahn-Hilliard models, Comput. Math. Appl., 68 (2014), 821–846. https://doi.org/10.1016/j.camwa.2014.07.014 doi: 10.1016/j.camwa.2014.07.014
![]() |
[2] |
J. Cahn, J. Hilliard, Free energy of a nonuniform system. Ⅰ. interfacial free energy, J. Chem. Phys., 28 (1958), 258–267. https://doi.org/10.1063/1.1744102 doi: 10.1063/1.1744102
![]() |
[3] |
A. Karma, W. Rappel, Quantitative phase-field modeling of dendritic growth in two and three dimensions, Phys. Rev. E, 57 (1998), 4323–4349. https://doi.org/10.1103/PhysRevE.57.4323 doi: 10.1103/PhysRevE.57.4323
![]() |
[4] |
S. Allen, J. Cahn, A microscopic theory for antiphase boundary motion and its application to antiphase domain coarsening, Acta. Matall., 27 (1979), 1085–1095. https://doi.org/10.1016/0001-6160(79)90196-2 doi: 10.1016/0001-6160(79)90196-2
![]() |
[5] |
R. Kobayashi, Modeling and numerical simulations of dendritic crystal growth, Physica D Nonlinear Phenom., 63 (1993), 410–423. https://doi.org/10.1016/0167-2789(93)90120-P doi: 10.1016/0167-2789(93)90120-P
![]() |
[6] |
M. Gurtin, D. Polignone, J. Vinals, Two-phase binary fluids and immiscible fluids described by an order parameter, Math. Models Methods Appl. Sci., 6 (1996), 815–831. https://doi.org/10.1142/S0218202596000341 doi: 10.1142/S0218202596000341
![]() |
[7] |
J. Barret, J. Blowey, H. Garcke, Finite element approximation of the Cahn-Hilliard equation with degenerate mobility, SIAM J. Numer. Anal., 37 (1999), 286–318. https://doi.org/10.1137/S0036142997331669 doi: 10.1137/S0036142997331669
![]() |
[8] |
C. Elliott, D. French, A nonconforming finite-element method for the two-dimensional Cahn-Hilliard equation, SIAM J. Numer. Anal., 26 (1989), 884–903. https://doi.org/10.1137/0726049 doi: 10.1137/0726049
![]() |
[9] |
C. Elliott, D. French, Numerical studies of the Cahn-Hilliard equation for phase separation, IMA J. Appl. Math., 38 (1987), 97–128. https://doi.org/10.1093/imamat/38.2.97 doi: 10.1093/imamat/38.2.97
![]() |
[10] |
J. Shen, J. Xu, J. Yang, The scalar auxiliary variable (SAV) approach for gradient flows, J. Comput. Phys., 353 (2018), 407–416. https://doi.org/10.1016/j.jcp.2017.10.021 doi: 10.1016/j.jcp.2017.10.021
![]() |
[11] |
S. Zhao, X. Xiao, X. Feng, An efficient time adaptivity based on chemical potential for surface Cahn-Hilliard equation using finite element approximation, Appl. Math. Comput., 369 (2020), 124901. https://doi.org/10.1016/j.amc.2019.124901 doi: 10.1016/j.amc.2019.124901
![]() |
[12] |
J. Shen, X. Yang, Numerical approximations of Allen-Cahn and Cahn-Hilliard equations, Discrete Contin. Dyn. Syst., 28 (2010), 1669–1691. https://doi.org/10.3934/dcds.2010.28.1669 doi: 10.3934/dcds.2010.28.1669
![]() |
[13] |
Y. Huang, W. Yang, H. Wang, J. Cui, Adaptive operator splitting finite element method for Allen-Cahn equation, Numer. Methods Partial Differ. Equations, 35 (2019), 1290–1300. https://doi.org/10.1002/num.22350 doi: 10.1002/num.22350
![]() |
[14] |
D. Kay, A. Tomasi, Color image segmentation by the Vector-Valued Allen-Cahn Phase-Field Model: A multigrid solution, IEEE. Trans. Image Process., 18 (2009), 2330–2339. https://doi.org/10.1109/TIP.2009.2026678 doi: 10.1109/TIP.2009.2026678
![]() |
[15] |
X. Feng, Y. Li, Analysis of interior penalty discontinuous Galerkin methods for the Allen-Cahn equation and the mean curvature flow, IMA J. Numer. Anal., 35 (2015), 1622–1651. https://doi.org/10.1093/imanum/dru058 doi: 10.1093/imanum/dru058
![]() |
[16] |
Y. Chen, Y. Huang, N. Yi, A SCR-based error estimation and adaptive finite element method for the Allen-Cahn equation, Comput. Math. Appl., 78 (2019), 204–223. https://doi.org/10.1016/j.camwa.2019.02.022 doi: 10.1016/j.camwa.2019.02.022
![]() |
[17] |
Y. Chen, Y. Huang, N. Yi, A decoupled energy stable adaptive finite element method for Cahn-Hilliard-Navier-Stokes equations, Commun. Comput. Phys., 29 (2021), 1186–1212. https://doi.org/10.4208/cicp.OA-2020-0032 doi: 10.4208/cicp.OA-2020-0032
![]() |
[18] |
D. Mao, L. Shen, A. Zhou, Adaptive finite element algorithms for eigenvalue problems based on local averaging type a posteriori error estimates, Adv. Comput. Math., 25 (2006), 135–160. https://doi.org/10.1007/s10444-004-7617-0 doi: 10.1007/s10444-004-7617-0
![]() |
[19] |
Z. Zhang, Z. Qiao, An adaptive time-stepping Strategy for the Cahn-Hilliard Equation, Commun. Comput. Phys., 11 (2012), 1261–1278. https://doi.org/10.4208/cicp.300810.140411s doi: 10.4208/cicp.300810.140411s
![]() |
[20] |
Y. Li, Y. Choi, J. Kim, Computationally efficient adaptive time step method for the Cahn-Hilliard equation, Comput. Math. Appl., 73 (2017), 1855–1864. https://doi.org/10.1016/j.camwa.2017.02.021 doi: 10.1016/j.camwa.2017.02.021
![]() |
[21] |
Y. Huang, N. Yi, The superconvergent cluster recovery method, J. Sci. Comput., 44 (2010), 301–322. https://doi.org/10.1007/s10915-010-9379-9 doi: 10.1007/s10915-010-9379-9
![]() |
[22] | A. Diegel, C. Wang, S. Wise, Stability and convergence of a second order mixed finite element method for the Cahn-Hilliard equation, arXiv preprint, 2016, arXiv: 1411.5248. https://doi.org/10.48550/arXiv.1411.5248 |
[23] |
J. Guo, C. Wang, S. Wise, X. Yue, An H2 convergence of a second-order convex-splitting, finite difference scheme for the three-dimensional Cahn–Hilliard equation, Commun. Math. Sci., 14 (2016), 489–515. https://dx.doi.org/10.4310/CMS.2016.v14.n2.a8 doi: 10.4310/CMS.2016.v14.n2.a8
![]() |
[24] |
K. Cheng, C. Wang, S. Wise, X. Yue, A second-order, weakly energy-stable pseudo-spectral scheme for the Cahn–Hilliard equation and its solution by the homogeneous linear iteration method, J. Sci. Comput., 69 (2016), 1083–1114. https://doi.org/10.1007/s10915-016-0228-3 doi: 10.1007/s10915-016-0228-3
![]() |
[25] |
J. Guo, C. Wang, S. Wise, X. Yue, An improved error analysis for a second-order numerical scheme for the Cahn-Hilliard equation, J. Comput. Appl. Math., 388 (2021), 113300. https://doi.org/10.1016/j.cam.2020.113300 doi: 10.1016/j.cam.2020.113300
![]() |
[26] | J. Shen, On error estimates of the projection methods for the Navier-Stokes equations: first-order schemes, Math. Comput., 65 (1996), 1039–1065. |
[27] |
X. Feng, A. Prohl, Error analysis of a mixed finite element method for the Cahn-Hilliard equation, Numer. Math., 99 (2004), 47–84. https://doi.org/10.1007/s00211-004-0546-5 doi: 10.1007/s00211-004-0546-5
![]() |
[28] |
C. Li, Y. Huang, N. Yi, An unconditionally energy stable second order finite element method for solving the Allen–Cahn equation, J. Comput. Appl. Math., 353 (2019), 38–48. https://doi.org/10.1016/j.cam.2018.12.024 doi: 10.1016/j.cam.2018.12.024
![]() |
[29] | Y. Yan, W. Chen, C. Wang, S. Wise, A second-order energy stable BDF numerical scheme for the Cahn-Hilliard equation, Commun. Comput. Phys., 23 (2018), 572–602. |
[30] |
K. Cheng, W. Feng, C. Wang, S. Wise, An energy stable fourth order finite difference scheme for the Cahn-Hilliard equation, J. Comput. Appl. Math., 362 (2019), 574–595. https://doi.org/10.1016/j.cam.2018.05.039 doi: 10.1016/j.cam.2018.05.039
![]() |
1. | Youngjin Hwang, Seokjun Ham, Chaeyoung Lee, Gyeonggyu Lee, Seungyoon Kang, Junseok Kim, A simple and efficient numerical method for the Allen–Cahn equation on effective symmetric triangular meshes, 2023, 31, 2688-1594, 4557, 10.3934/era.2023233 |
h | ‖umh−umh2‖0 | Rate | △t | Rate |
116 | 0.052470400 | 116 | ||
132 | 0.018185700 | 1.528700 | 132 | 2.18632 |
164 | 0.003449260 | 2.398450 | 164 | 2.52708 |
1128 | 0.000945365 | 1.867340 | 1128 | 1.78946 |