Loading [MathJax]/jax/output/SVG/jax.js
Special Issues

High-order energy stable schemes of incommensurate phase-field crystal model

  • This article focuses on the development of high-order energy stable schemes for the multi-length-scale incommensurate phase-field crystal model which is able to study the phase behavior of aperiodic structures. These high-order schemes based on the scalar auxiliary variable (SAV) and spectral deferred correction (SDC) approaches are suitable for the L2 gradient flow equation, i.e., the Allen-Cahn dynamic equation. Concretely, we propose a second-order Crank-Nicolson (CN) scheme of the SAV system, prove the energy dissipation law, and give the error estimate in the almost periodic function sense. Moreover, we use the SDC method to improve the computational accuracy of the SAV/CN scheme. Numerical results demonstrate the advantages of high-order numerical methods in numerical computations and show the influence of length-scales on the formation of ordered structures.

    Citation: Kai Jiang, Wei Si. High-order energy stable schemes of incommensurate phase-field crystal model[J]. Electronic Research Archive, 2020, 28(2): 1077-1093. doi: 10.3934/era.2020059

    Related Papers:

    [1] Kai Jiang, Wei Si . High-order energy stable schemes of incommensurate phase-field crystal model. Electronic Research Archive, 2020, 28(2): 1077-1093. doi: 10.3934/era.2020059
    [2] 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
    [3] 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
    [4] Junseok Kim . Maximum principle preserving the unconditionally stable method for the Allen–Cahn equation with a high-order potential. Electronic Research Archive, 2025, 33(1): 433-446. doi: 10.3934/era.2025021
    [5] Xiang Xu . Recent analytic development of the dynamic $ Q $-tensor theory for nematic liquid crystals. Electronic Research Archive, 2022, 30(6): 2220-2246. doi: 10.3934/era.2022113
    [6] Wenya Qi, Padmanabhan Seshaiyer, Junping Wang . A four-field mixed finite element method for Biot's consolidation problems. Electronic Research Archive, 2021, 29(3): 2517-2532. doi: 10.3934/era.2020127
    [7] Cui Guo, Yixue Wang, Haibin Wang, Xiongbo Zheng, Bin Zhao . Non-destructive test method of wood moisture content based on multiple varying bounds integral numerical method. Electronic Research Archive, 2025, 33(4): 2246-2274. doi: 10.3934/era.2025098
    [8] Cheng Wang . Convergence analysis of Fourier pseudo-spectral schemes for three-dimensional incompressible Navier-Stokes equations. Electronic Research Archive, 2021, 29(5): 2915-2944. doi: 10.3934/era.2021019
    [9] Li Tian, Ziqiang Wang, Junying Cao . A high-order numerical scheme for right Caputo fractional differential equations with uniform accuracy. Electronic Research Archive, 2022, 30(10): 3825-3854. doi: 10.3934/era.2022195
    [10] Leilei Wei, Xiaojing Wei, Bo Tang . Numerical analysis of variable-order fractional KdV-Burgers-Kuramoto equation. Electronic Research Archive, 2022, 30(4): 1263-1281. doi: 10.3934/era.2022066
  • This article focuses on the development of high-order energy stable schemes for the multi-length-scale incommensurate phase-field crystal model which is able to study the phase behavior of aperiodic structures. These high-order schemes based on the scalar auxiliary variable (SAV) and spectral deferred correction (SDC) approaches are suitable for the L2 gradient flow equation, i.e., the Allen-Cahn dynamic equation. Concretely, we propose a second-order Crank-Nicolson (CN) scheme of the SAV system, prove the energy dissipation law, and give the error estimate in the almost periodic function sense. Moreover, we use the SDC method to improve the computational accuracy of the SAV/CN scheme. Numerical results demonstrate the advantages of high-order numerical methods in numerical computations and show the influence of length-scales on the formation of ordered structures.



    Aperiodic crystals, such as quasicrystals, are an important class of materials whose Fourier spectra cannot be all expressed by a set of basis vectors over the rational number field. The irrational coefficients give rise to the denseness of Fourier spectra which results in the difficulties in the theoretical study. Theoretically, a multiple characteristic length-scale model which possesses, at least, an irrational scale, has been widely applied to study the formation and thermodynamic stability of the aperiodic structures [1, 6, 13, 10, 14, 7]. The early model could trace back to Bak's work on three-dimensional icosahedral quasicrystals. Since then, many related models have been proposed to study aperiodic structures, including for multicomponent systems [7]. Among these models, Lifshitz and Petrich (LP) modified the Swift-Hohenberg model and explicitly added an incommensurate two-length-scale potential into a Lyapunov functional to explore quasiperiodic patterns that emerged in Faraday experiments [13]. Recently, Savitz et al. extended the LP model from two-length-scale potential to multiple (3) length-scale potential and studied more kinds of quasicrystals [14]. The m-length-scale model actually is an incommensurate multi-length-scale phase-field-crystal (iPFC) model who owns a 4m (mN) order differential operator and nonlinear term in the energy functional. A high precision computation is helpful to study the phase behaviors of aperiodic structures. In this article, we will pay attention to the development of high-order numerical methods for the iPFC model.

    Recently, various numerical methods have been proposed to solve phase-field equations including the convex splitting methods [17], the linear stabilized schemes [16], the invariant energy quadratization (IEQ) [18] and the scalar auxiliary variable (SAV) approaches [15,11]. The convex splitting method splits the energy functional into the convex and concave parts. The method treats the convex part implicitly and the concave one explicitly to keep the unconditional energy stability. While the application of this method is restricted by the form of the energy functional, such as double-well bulk energy. The linear stabilized scheme adds a penalty term to improve its stability and deals with the nonlinear terms explicitly for implementing it easily. However, such a stabilized approach makes it difficult to design second-order unconditionally energy stable schemes. Assume that the nonlinear part has a lower bound, via introducing an auxiliary variable, the IEQ method transforms the energy functional into a quadratic form to keep the unconditional energy dissipation property. Similarly, the SAV approach introduces a scalar auxiliary variable by supposing the bounded bulk energy and obtains an unconditionally energy stable system. Besides these methods, the spectral deferred correction (SDC) [4,5] algorithm is an efficient strategy to improve the accuracy of the above schemes. In the paper, we will apply the SAV approach to solve the time-dependent equations and further use the SDC strategy to improve the numerical accuracy.

    For aperiodic structures, two kinds of numerical methods, including the crystalline approximant method (CAM) and the projection method (PM) are usually used to discretize the quasiperiodic functions [9]. The CAM uses a big periodic structure to approximate an aperiodic structure and corresponds to the Diophantine approximation problem which studies how to approximate irrational numbers by rational numbers [3]. To evaluate aperiodic structures accurately, the CAM needs an extremely big computational region with an unacceptable computational burden to reduce the error of Diophantine approximation. To avoid the Diophantine approximation problem, the PM accurately describes aperiodic structures based on the fact that the aperiodic structure can be regarded as a periodic crystal in an appropriate higher-dimensional space. The PM uses one higher-dimensional periodic region to capture the essential characteristics and greatly reduces computational complexity.

    The rest of the paper is organized as follows. In Section 2, we first outline some useful preliminaries of the almost periodic functions and then present the iPFC model. The energy dissipation principle of the L2 gradient flow for the iPFC model in the almost periodic sense is also given. In Section 3, we propose a second-order energy stable scheme for the iPFC model and give the corresponding error estimate. And we use the SDC approach to improve its computational accuracy efficiently. Section 4 presents the convergence rates of these numerical schemes and discusses the advantages of high-order numerical approaches in simulating dynamic evolution. Moreover, we also show the influence of multiple length-scales on the thermodynamic stability of aperiodic structures. There are some conclusions in Section 5.

    Aperiodic structures are space-filling phases without decay. A useful mathematical theory to describe aperiodic structures is the almost periodic function theory which is a generalization of continuous periodic functions. We define the notation of a d-dimensional almost periodic function.

    Definition 2.1. Let f(r) be a real-valued or complex-valued function defined on Rd and let ϵ>0. We say that ζRd is an ϵ-almost period of f if

    |f(rζ)f(r)|<ϵ, for   any  rRd.

    A function f is almost periodic on Rd if it is continuous and if for every ϵ there exists a number L=L(ϵ,f) such that any cube with the side length of L on Rd contains an ϵ-almost period of f.

    The almost periodic L2 inner product is defined by

    f,gAP=limR1|Q(R)|Q(R)f(r)¯g(r)dr,

    where Q(R)=[R,R]d and |Q(R)| is the measure of Q(R). It is known that this inner product is well-defined [2]. We denote the corresponding norm by 2AP=,AP. For simplicity, we denote the average spacial integral over the whole space as

    =limR1|Q(R)|Q(R).

    Some useful properties of d-dimensional almost periodic functions are presented as follows.

    Proposition 1. (1). An almost periodic function is uniformly continuous and bounded.

    (2). If f(r) and g(r) are almost periodic functions, then f(r)+g(r) and f(r)g(r) are almost periodic functions.

    These properties can be easily proven from the one-dimensional results [2].

    Theorem 2.2. If f(r) and g(r) are almost periodic functions and differentiable, then the Green's identity holds in the almost periodic sense, i.e.,

    f(r),g(r)AP=f(r),g(r)AP.

    Proof. Since f(r) and g(r) are almost periodic functions, then there exists a constant M such that supr{|f(r)|,|g(r)|}M. Denote

    bR=limR1|Q(R)|Q(R)f(r)n¯g(r)ds,

    where n is the outward normal of Q(R). In the d-dimensional space, we have

    |bR|limRM22d(2R)d1(2R)d=0.

    Therefore, we obtain the desired conclusion by

    f(r),g(r)AP=f(r)g(r)dr=f(r)g(r)dr=f(r),g(r)AP.

    The simplest iPFC model may be the LP model which was originally proposed to study the bi-frequency excited Faraday wave [13]. Concretely, the free energy functional of the LP model can be written as

    FLP[ψ(r)]={c2[(Δ+1)(Δ+q2)ψ]2+(ε2ψ2α3ψ3+14ψ4)}dr,

    where q is an irrational number depending on the property of quasiperiodic structures. The essential feature of the energy functional is the existence of two characteristic length-scales, 1 and q, which is a critical factor to stabilize quasi-periodic structures. The LP model has been also used to study the soft-matter quasicrystals [12]. However, this model could be able to globally stabilize the dodecagonal and decagonal rotational quasicrystals [8]. To generalize the iPFC model, Savitz et al. extended the interaction potential from two characteristic length-scales to multiple characteristic length-scales. In particular, the Lyapunov functional of an m-length-scale incommensurate system can be written as

    F[ψ(r)]={c2[mj=1(Δ+q2j)ψ]2+(ε2ψ2α3ψ3+14ψ4)}dr,

    where ψ(r), rRd (d=1,2,3), is the order parameter corresponding to the density profile of the system. qj is the j-th characteristic length-scale which depends on the property of aperiodic structures. c is a positive model parameter to ensure that the principle wavelengths are near the critical wavelengths. ε and α are both model parameters related to physical conditions, such as temperature, pressure. To reduce the number of model parameters, the iPFC energy functional can be rescaled by defining F=c2F, ψ(r)=cϕ(r), ˜ε=cε, and ˜α=cα. The rescaled functional is

    F[ϕ(r)]={12[Gϕ]2+N(ϕ)}dr=12Gϕ2AP+N(ϕ),1AP,

    where

    G=mj=1(Δ+q2j),N(ϕ)=˜ε2ϕ2˜α3ϕ3+14ϕ4.

    To solve the iPFC free energy functional, we consider the following Allen-Cahn dynamic equation

    ϕt=W(ϕ),W(ϕ):=δFδϕ=G2ϕ+N(ϕ). (1)

    The initial value is ϕ|t=0=ϕ0. From Theorem 2.2, we can prove that the system (11) satisfies the following energy dissipation law

    dF(ϕ)dt=δFδϕ,ϕtAP=W(ϕ)2AP0.

    Then we impose the following mean zero constraint of order parameter on the iPFC model to ensure the mass conservation

    ϕ(r)dr=0. (2)

    In this section, we propose the second-order unconditionally energy stable Crank-Nicolson scheme based on the SAV technique, and also give the error analysis. Further, we use the SDC strategy to improve the temporal accuracy of the second-order scheme.

    Let's suppose that F1(ϕ)=N(ϕ),1AP+C10, where C1 is a positive constant. We introduce a scalar auxiliary variable R=F1(ϕ) to transform (1) into an equivalent system as

    ϕt=W(ϕ), (3a)
    W(ϕ)=G2ϕ+RF1(ϕ)N(ϕ), (3b)
    Rt=N(ϕ)2F1(ϕ),ϕtAP. (3c)

    By taking the almost periodic inner products of (3a) with W, and (3b) with ϕt, multiplying (3c) with 2R and adding them together, we have the following energy dissipation property

    ddt(12Gϕ2AP+R2C1)=W2AP0.

    The SAV approach can construct high-order unconditionally energy stable schemes. In this section, we discuss a second-order semi-discrete scheme based on the Crank-Nicolson method. Suppose the time interval [0,T] is divided into NT non-overlapping subintervals by the partition 0=t0<t1<<tn<<tNT=T. The time step size is τn=tn+1tn. We denote tn+τn/2 as tn+1/2. Let ϕn=ϕ(tn) and Rn=R(tn).

    Scheme 3.1 (SAV/CN). For n1, given ϕn, Rn and ϕn1, Rn1, we update ϕn+1 and Rn+1 by

    ϕn+1ϕnτn=Wn+1/2, (4a)
    Wn+1/2=G2ϕn+1/2+Rn+1/2F1(ˉϕn+1/2)N(ˉϕn+1/2), (4b)
    Rn+1Rnτn=N(ˉϕn+1/2)2F1(ˉϕn+1/2),ϕn+1ϕnτnAP, (4c)

    where Rn+1/2=(Rn+1+Rn)/2, ϕn+1/2=(ϕn+1+ϕn)/2 and ˉϕn+1/2=(3ϕnϕn1)/2.

    Theorem 3.1. The SAV/CN scheme satisfies the following energy dissipation mechanism

    Fn+1SAV/CNFnSAV/CN0,

    where

    FnSAV/CN=12Gϕn2AP+(Rn)2C1.

    Proof. We take the almost periodic inner products of (4a) with τnWn+1/2, and (4b) with (ϕn+1ϕn).

    ϕn+1ϕn,Wn+1/2AP=Wn+1/22AP, (5)
    ϕn+1ϕn,Wn+1/2AP=ϕn+1ϕn,G2ϕn+1/2APϕn+1ϕn,Rn+1/2F1(ˉϕn+1/2)N(ˉϕn+1/2)AP. (6)

    Multiplying (4c) with 2Rn+1/2 yields

    2Rn+1/2(Rn+1Rn)=Rn+1/2F1(ˉϕn+1/2)N(ˉϕn+1/2),ϕn+1ϕnAP. (7)

    Adding (5), (6) and (7) together, we obtain

    2Rn+1/2(Rn+1Rn)=ϕn+1ϕn,G2ϕn+1/2APWn+1/22AP. (8)

    Since Rn+1/2=(Rn+1+Rn)/2 and ϕn+1/2=(ϕn+1+ϕn)/2, (8) can be recast as

    12(Gϕn+12APGϕn2AP)+(Rn+1)2(Rn)2=Wn+1/22AP0.

    The desired conclusion is obtained from the above equation.

    Remark 3.1. It is noted that the modified energy FnSAV/CN is different from the original energy F(ϕn) since Rn is obtained from the iteration process.

    Remark 3.2(The implementation of the SAV/CN scheme). Denote

    u(tn+1/2)=N(ϕ(tn+1/2))F1(ϕ(tn+1/2)),un+1/2=N(ˉϕn+1/2)F1(ˉϕn+1/2).

    Substituting (4b) and (4c) into (4a), we obtain

    ϕn+1ϕnτn=[G2ϕn+1/2+un+1/2(Rn+14un+1/2,ϕn+1ϕnAP)]. (9)

    Eqn. (9) can be rewritten as

    (I+τn2G2)ϕn+1+τn4un+1/2un+1/2,ϕn+1AP=(Iτn2G2)ϕnτnRnun+1/2+τn4un+1/2,ϕnAPun+1/2. (10)

    Taking the almost periodic inner product with (I+12τnG2)1un+1/2 leads to

    un+1/2,ϕn+1AP+τn4γnun+1/2,ϕn+1AP=un+1/2,(I+τn2G2)1cnAP, (11)

    where

    γn=un+1/2,(I+τn2G2)1un+1/2AP,
    cn=(Iτn2G2)ϕnτnRnun+1/2+τn4un+1/2,ϕnAPun+1/2.

    From (11), we get

    un+1/2,ϕn+1AP=un+1/2,(I+12τnG2)1cnAPI+14τnγn. (12)

    Then we can directly calculate ϕn+1 via (10) and (12).

    In this section, we will derive the error estimate of the SAV/CN scheme 3.1. Denote en=ϕnϕ(tn), wn+1=Wn+1W(tn+1), and rn=RnR(tn), we have

    Theorem 3.2. For the Allen-Cahn dynamic equation, assume that ϕ0 is almost periodic and ϕt2AP is bounded. Considering a uniform time partition, i.e., τ=τk,kNT, we have

    Gek2AP+(rk)2Cτ4tk0(ϕttt(s)2AP+|rttt(s)|2)ds,

    where the constant C is independent on τ.

    Proof. Let's subtract (3) from (4) at tn+1/2

    en+1en=τwn+1/2+Tn+1/21, (13)
    wn+1/2=G2en+1/2+Rn+1/2un+1/2R(tn+1/2)u(tn+1/2), (14)
    rn+1rn=12un+1/2,ϕn+1ϕnAP12u(tn+1/2),τϕt(tn+1/2)AP+Tn+1/22, (15)

    where the truncation errors are given by

    Tn+1/21=τϕt(tn+1/2)(ϕ(tn+1)ϕ(tn)),
    Tn+1/22=τRt(tn+1/2)(R(tn+1)R(tn)).

    With the Taylor expansion, the truncation errors can be rewritten as

    Tn+1/21=12tn+1/2tn+1(tn+1s)2ϕttt(s)ds12tn+1/2tn(tns)2ϕttt(s)ds,
    Tn+1/22=12tn+1/2tn+1(tn+1s)2Rttt(s)ds12tn+1/2tn(tns)2Rttt(s)ds.

    Firstly, making the almost periodic inner product of (13) with wn+1/2 yields

    en+1en,wn+1/2AP+τwn+1/22AP=Tn+1/21,wn+1/2AP. (16)

    Then its right-hand term can be bounded by

    Tn+1/21,wn+1/2APτ2wn+1/22AP+CτTn+1/212APτ2wn+1/22AP+Cτ4tn+1tnϕttt(s)2APds.

    Secondly, by taking the almost periodic inner products of (14) with (en+1en), we obtain

    wn+1/2,en+1enAP=12(Gen+12APGen2AP)Rn+1/2un+1/2R(tn+1/2)u(tn+1/2),en+1enAP. (17)

    Without the minus sign, the second term on the right-hand side in the above equation can be transformed into

    Rn+1/2un+1/2R(tn+1/2)u(tn+1/2),en+1enAP=rn+1/2un+1/2,en+1enAP+R(tn+1/2)un+1/2u(tn+1/2),en+1enAP. (18)

    Note that |R(t)|C and

    un+1/2u(tn+1/2)APCN(ˉϕn+1/2)N(ϕ(tn+1/2))APC(enAP+en1AP).

    The last term on the right-hand side of (18) can be estimated by

    R(tn+1/2)un+1/2u(tn+1/2),en+1enAP=R(tn+1/2)un+1/2u(tn+1/2),τwn+1/2+Tn+1/21APτ2wn+1/22AP+Cτun+1/2u(tn+1/2)2AP+CτTn+1/212APτ2wn+1/22AP+Cτ(en2AP+en12AP)+Cτ4tn+1tnϕttt(s)2APds.

    Thirdly, multiplying the both sides of (15) by 2rn+1/2, we obtain

    (rn+1)2(rn)2=rn+1/2un+1/2,ϕn+1ϕnAPrn+1/2u(tn+1/2),τϕt(tn+1/2)AP+2rn+1/2Tn+1/22. (19)

    The first two terms on the right-hand side of (19) can be rewritten as

    rn+1/2un+1/2,ϕn+1ϕnAPrn+1/2u(tn+1/2),τϕt(tn+1/2)AP=rn+1/2un+1/2,en+1enAPrn+1/2u(tn+1/2),Tn+1/21AP+rn+1/2un+1/2u(tn+1/2),ϕ(tn+1)ϕ(tn)AP,

    where the last two terms on the right-hand side satisfy

    rn+1/2u(tn+1/2),Tn+1/21APCτ((rn+1)2+(rn)2)+Cτ4tn+1tnϕttt(s)2APds,

    and

    rn+1/2un+1/2u(tn+1/2),ϕ(tn+1)ϕ(tn)APCτϕt2AP((rn+1/2)2+un+1/2u(tn+1/2)2AP)Cτ((rn+1)2+(rn)2+en2AP+en12AP).

    Therefore the last term on the right-hand side of (19) can be bounded by

    2rn+1/2Tn+1/22Cτ((rn+1)2+(rn)2)+Cτ4tn+1tn|rttt(s)|2ds.

    With the above estimates, adding (16), (17) and (19) together leads to

    12(Gen+12APGen2AP)+(rn+1)2(rn)2Cτ(en2AP+en12AP+(rn+1)2+(rn)2)+Cτ4tn+1tnϕttt(s)2APds+Cτ4tn+1tn|rttt(s)|2ds.

    Then we obtain the desired conclusion by summing over n, n=0,1,,k1, and using the Gronwall inequality.

    The SDC approach [4,5] is an efficient method to improve the numerical precision for an existing scheme using spectral collocation points. Firstly, we introduce the basic idea of the SDC method. Integrating both sides of the Eqn.(1) with respect to t, we have

    ϕ(t)=ϕ(0)t0W(ϕ(τ))dτ. (20)

    Suppose the approximation solution ϕ[0](t) has been calculated by some numerical schemes. Then we define the residual R[0](t) and the error ϵ[0](t) as

    R[0](t)=ϕ(0)t0W(ϕ[0](τ))dτϕ[0](t), (21)
    ϵ[0](t)=ϕ(t)ϕ[0](t). (22)

    Replacing (22) into (20) yields

    ϕ(t)=ϕ(0)t0W(ϕ[0](τ)+ϵ[0](τ))dτ. (23)

    Then we insert (23) into (22) and subtract (21)

    ϵ[0](t)R[0](t)=t0W(ϕ[0](τ)+ϵ[0](τ))dτ+t0W(ϕ[0](τ))dτ.

    By taking the derivative of both sides of the above equation, we obtain

    dϵ[0](t)dt=W(ϕ[0](t)+ϵ[0](t))+W(ϕ[0](t))+dR[0](t)dt. (24)

    Then, we apply the SDC approach into the second-order SAV/CN scheme to improve the numerical accuracy. We denote this strategy as SAV/CN+SDC. To calculate the integral precisely, we adopt the following Chebyshev nodes,

    tn=T2T2cos(nπNT),n=0,1,,NT.

    The time step size is τn=tn+1tn. To obtain the approximation solution ϕ[0](t), we rewrite the SAV/CN scheme as

    ϕn+1[0]ϕn[0]τn=Wn+1/2[0],Wn+1/2[0]=G2ϕn+1/2[0]+Rn+1/2[0]F1(ˉϕn+1/2[0])N(ˉϕn+1/2[0]),Rn+1[0]Rn[0]τn=N(ˉϕn+1/2[0])2F1(ˉϕn+1/2[0]),ϕn+1[0]ϕn[0]τnAP, (25)

    where ϕ0[0]=ϕ0, Rn+1/2[0]=(Rn+1[0]+Rn[0])/2, ϕn+1/2[0]=(ϕn+1[0]+ϕn[0])/2 and ˉϕn+1/2[0]=(3ϕn[0]ϕn1[0])/2. We adopt a similar strategy to discretize (24) as follows

    ϵn+1[0]ϵn[0]τn=Wn+1/2,ϵ[0]+Wn+1/2[0]+Rn+1[0]Rn[0]τn, (26a)
    Wn+1/2,ϵ[0]=G2ϕn+1/2,ϵ[0]+Rn+1/2[0]F1(ˉϕn+1/2[0])N(ˉϕn+1/2,ϵ[0]), (26b)

    where ϵ0[0]=ϕ(0)ϕ0[0]=0, ϕn+1/2,ϵ[0]=ϕn+1/2[0]+ϵn+1/2[0], and ˉϕn+1/2,ϵ[0]=ˉϕn+1/2[0]+ˉϵn+1/2[0]. Substituting (25) and (26b) into (26a) and eliminating the residual by (21), we obtain the following linear solvable equation

    (I+τn2G2)ϵn+1[0]=(Iτn2G2)ϵn[0]tn+1tnW(ϕ[0](τn))dτnϕn+1[0]+ϕn[0]τnRn+1/2[0]F1(ˉϕn+1/2[0]){N(ˉϕn+1/2,ϵ[0])N(ˉϕn+1/2[0])}.

    A more accurate solution can be updated by

    ϕn+1[1]=ϕn+1[0]+ϵn+1[0].

    The PM is an accurate approach in computing aperiodic structures that can avoid the Diophantine approximation error. The PM is based on the fact that the d-dimensional aperiodic structure can be embedded into an n-dimensional periodic structure (nd). The dimensionality n is determined by the spectrum structures of the aperiodic structures. In particular, n is the number of linearly independent vectors over the rational number field which span the spectrum. In the PM, the d-dimensional order parameter ϕ(r) can be expanded as

    ϕ(r)=hZnˆϕ(h)ei[(PBh)Tr],rRd.

    where BRn×n is an invertible matrix related to the n-dimensional primitive reciprocal lattice. The projection matrix PRd×n depends on the property of aperiodic structures. If we consider d-dimensional periodic structures, the projection matrix degenerates to a d-order identity matrix. Therefore, the PM provides a unified framework in calculating periodic and aperiodic crystals. More details about the PM can refer to [9]. The Fourier coefficient ˆϕ(h) satisfies

    X:={(ˆϕ(h))hZn:ˆϕ(h)C,hZn|ˆϕ(h)|<}.

    In practice, let N=(N1,N2,,Nn)Nn, and

    XN:={ˆϕ(h)X:ˆϕ(h)=0,forall|hj|>Nj/2,j=1,2,,n}.

    The number of elements in the set is N=(N1+1)(N2+1)(Nn+1). Using the PM, the SAV/CN scheme of full discretization reads

    ˆϕn+1(h)ˆϕn(h)=τˆWnh(h),ˆWnh(h)=12mj=1[q2j(PBh)T(PBh)]2[ˆϕn+1(h)+ˆϕn(h)]+Rn+1+Rn2Fnh1[ˆΦ]^Nnh(h),Rn+1Rn=h1+h2=0^Nnh(h1)2Fnh1[ˆΦ][ˆϕn+1(h2)ˆϕn(h2)],

    where

    ^Nnh(h)=˜εˆϕnh(h)˜αh1+h2=hˆϕnh(h1)ˆϕnh(h2)+h1+h2+h3=hˆϕnh(h1)ˆϕnh(h2)ˆϕnh(h3),Fnh1[ˆΦ]=˜ε2h1+h2=0ˆϕnh(h1)ˆϕnh(h2)˜α3h1+h2+h3=0ˆϕnh(h1)ˆϕnh(h2)ˆϕnh(h3)+14h1+h2+h3+h4=0ˆϕnh(h1)ˆϕnh(h2)ˆϕnh(h3)ˆϕnh(h4)+C1,ˆϕnh(h)=3ˆϕn+1(h)ˆϕn(h)2.

    In the above equations, the nonlinear terms are n-dimensional convolutions in the Fourier space. Directly computing them is extremely expensive. To avoid this, we use the pseudospectral method through the n-dimensional fast Fourier transformation to compute them efficiently in the n-dimensional time domain by simple multiplication. The mass conservation constraint (2) can be satisfied through

    eT1ˆΦ=0,

    where e1=(1,0,,0)TRN.

    In this section, we present several numerical examples to verify the accuracy of the SAV/CN and SAV/CN+SDC schemes and to illustrate the advantages of the higher-order scheme in dynamic evolution. We also show the influence of multiple length-scales on the thermodynamic stability of aperiodic structures.

    In this subsection, we take the two characteristic length scale iPFC model in one-dimensional space to test the numerical accuracy of the SAV/CN and SAV/CN+SDC schemes. The model parameters are set as q1=2, q2=3, ˜ε=10 and ˜α=4. The computational domain is [0,2π]. Correspondingly, the projection matrix P and B in the PM both are 1. The initial data is chosen as ϕ(x,0)=sin(x). We check the temporal accuracy by taking the space discretization N=128. The numerical solution of NT=2048 is set as the reference solution. Tab. 1 shows the errors and convergence rates of the SAV/CN and SAV/CN+SDC schemes at T=0.2. One can observe that the numerical accuracy of the SAV/CN scheme is second-order and can be improved to fourth-order by the SDC approach.

    Table 1.  Errors and convergence rates of the SAV/CN and SAV/CN+SDC schemes for the Allen-Cahn equation. The numerical solution of NT=2048 is regarded as the reference value.
    64 128 256 512
    SAV/CN Error 4.75E-3 1.17E-3 2.91E-4 7.17E-5
    Rate - 2.01 2.02 2.07
    SAV/CN + SDC Error 1.16E-5 6.78E-7 4.04E-8 2.46E-9
    Rate - 4.07 4.03 4.02

     | Show Table
    DownLoad: CSV

    In this subsection, we simulate the dynamic process of 2-dimensional dodecagonal quasicrystal (DDQC) using the iPFC model with two length-scales. The model parameters are q1=1, q2=2cos(π/12), ˜ε=2 and ˜α=2. The initial value is the DDQC whose spectral distribution and real morphology are shown in Fig. 1. The big blue dot represents the origin and the others are the 24 basic Fourier modes located on the circles of radii q1 and q2, respectively. When calculating the DDQC, we adopt the PM to discretize the spatial functions in 4-dimensional space with 244 trigonometric functions. The 244 basis functions bring a negligible spatial error comparing the temporal error. The projection matrix in the PM is

    P=(1cos(π/6)cos(π/3)00sin(π/6)sin(π/3)1),
    Figure 1.  The spectral distribution and the corresponding real morphology of the initial value.

    and the B is a 4-order identity matrix. We use the second-order SAV/CN scheme with NT=256 to simulate the dynamic evolution for the DDQC. Fig. 2 shows the change tendency of the energy value. The corresponding morphologies at t=0,50,100,150,200 are presented in Fig. 3. As one can see, the proposed scheme satisfies the energy dissipation law. It should be noted that the value C1 in the SAV approach plays an important role in computational simulations. In our computation, the C1 chosen as 1016 which maintains the consistency of original and modified energy values.

    Figure 2.  The time evolution of energy which is obtained by the SAV/CN scheme in the case of NT=256. The model parameters are q1=1, q2=2cos(π/12), ˜ε=2 and ˜α=2.
    Figure 3.  The morphologies of dynamic evolution in Fig. 2. Snapshots are taken at t=50,100,150,200, respectively.

    To show the role of the high-order methods in dynamic evolution, we give the reference energy Es which is calculated by the SAV/CN+SDC scheme in the case of NT=2048. Using the reference value as the baseline, Fig. 4 shows the energy difference of the SAV/CN scheme with NT=64,128,256 and SAV/CN+SDC method with NT=32. With the increase of the temporal discretization NT, the energy difference decreases. However, the energy value obtained by the fourth-order SAV/CN+SDC scheme with NT=32 is closer to the reference value than that of the SAV/CN scheme. It is demonstrated that the higher-order method shows more accurate results with less time discretization points in the dynamic simulation. We also give the morphologies of the crucial moment t=0.4025 in Fig. 5. The morphology of the reference solution is also set as a reference value to clearly illustrate the differences. And the conclusions from these results are consistent with the energy evolution curves.

    Figure 4.  The difference between the numerical energy values and the reference value Es. The numerical energy values are computed by the two methods: SAV/CN and SAV/CN+SDC. The reference energy value is obtained by the SAV/CN+SDC method in the case of NT=2048. The model parameters are q1=1, q2=2cos(π/12), ˜ε=2 and ˜α=2.

    In this subsection, we use the dodecagonal quasiperiodic phases as an example to investigate the influence of multi-length-scale potentials on the stability of aperiodic structures. Concretely, we consider three, four, and five multiple length scale iPFC models. The parameters in the potential function are set as qj=sj1, s=2cos(π/12), j=1,,m, m=3,4,5. The other model parameters are ˜ε=2 and ˜α=2. In the PM, the P and B are consistent with Subsec. 4.2. The spatial functions are also discretized in 4-dimensional space with 244 basis functions. The SAV/CN approach with NT=256 is adopted to simulate the dynamic process. Fig. 6 lists the corresponding energy evolution plots, initial values, and stationary states. The first row presents the energy evolutions of DDQCs with different length-scale potentials. The second and third rows give the spectral distributions and real morphologies of the m-length-scale DDQCs, respectively. The last row shows the real morphologies of stationary solutions. From these results, one can see that the three- and four-length-scale potentials both result in 6-fold symmetric periodic crystal, while the five-length-scale potential can obtain the 12-fold symmetric quasicrystals. Therefore increasing the number of characteristic length-scales in the iPFC model is helpful to stabilize the quasicrystals.

    Figure 5.  The real morphologies which describe the difference between the numerical solutions and the reference value at t=0.4025. The reference solution is obtained by the SDC/CN+SDC method in the case of NT=2048. The model parameters are set as q1=1, q2=2cos(π/12), ˜ε=2 and ˜α=2.
    Figure 6.  The energy evolutions, initial values and convergence solutions of m-length-scale DDQCs under the model parameters ˜ε=2, ˜α=2. The scale parameters are set as qj=sj1, j=1,,m, s=2cos(π/12).

    For the iPFC model, we proposed a second-order SAV/CN scheme which is unconditionally energy stable in the almost periodic function sense and gave the error estimate. Meanwhile, we used the SDC approach to further improve the accuracy of the second-order scheme to the fourth-order method through a one-step correction. The PM was applied to discretize spatial functions for computing aperiodic structures to high accuracy. In numerical simulations, the efficiency of the numerical schemes was demonstrated via the numerical convergence rates and the comparison of the dynamic evolutions. By comparing the dynamic evolutions of the DDQCs with different length-scales, we found that increasing the number of characteristic length-scales in the iPFC model significantly is helpful to stabilize aperiodic structures.

    This work is supported by National Natural Science Foundation of China (11771368), Hunan Science Foundation of China (2018JJ2376), Youth Project (18B057) and Key Project (19A500) of Education Department of Hunan Province of China, and Hunan Provincial Innovation Foundation for Postgraduate (0943/431000232).



    [1] Phenomenological theory of icosahedral incommensurate ("quasiperiodic") order in Mn-Al alloys. Phys. Rev. Lett. (1985) 54: 1517-1519.
    [2] C. Corduneanu, Almost Periodic Functions, 2nd edition, Chelsea Publishing Company, New York, 1989.
    [3] Simultaneous Diophantine approximation. Duke Math. J. (1946) 13: 105-111.
    [4] Spectral deferred correction methods for ordinary differential equations. BIT Numer. Math. (2000) 40: 241-266.
    [5] Semi-implicit spectral deferred correction method based on the invariant energy quadratization approach for phase field problems. Commun. Comput. Phys. (2019) 26: 87-113.
    [6] Long-range icosahedral orientational order and quasicrystals.. Phys. Rev. Lett. (1985) 55: 607-610.
    [7] Stability of three-dimensional icosahedral quasicrystals in multi-component systems. Philos. Mag. (2020) 100: 84-109.
    [8] K. Jiang, J. Tong, P. Zhang and A.-C. Shi, Stability of two-dimensional soft quasicrystals in systems with two length scales, Phys. Rev. E, 92 (2015), 042159.
    [9] Numerical methods for quasicrystals. J. Comput. Phys. (2014) 256: 428-440.
    [10] K. Jiang, P. Zhang and A.-C. Shi, Stability of icosahedral quasicrystals in a simple model with two-length scales, J. Phys.: Condens. Matter, 29 (2017), 124003. doi: 10.1088/1361-648X/aa586b
    [11] X. Li and J. Shen, Stability and error estimates of the SAV Fourier-spectral method for the phase field crystal equation, Advances in Computational Mathematics, arXiv: 1907.07462, 46 (2020), Article number: 48. doi: 10.1007/s10444-020-09789-9
    [12] Soft quasicrystals–Why are they stable?. Philos. Mag. (2007) 87: 3021-3030.
    [13] Theoretical model for Faraday waves with multiple-frequency forcing. Phys. Rev. Lett. (1997) 79: 1261-1264.
    [14] Multiple-scale structures: From Faraday waves to soft-matter quasicrystals. IUCrJ (2018) 5: 247-268.
    [15] The scalar auxiliary variable (SAV) approach for gradient flows. J. Comput. Phys. (2018) 353: 407-416.
    [16] Numerical approximations of Allen-Cahn and Cahn-Hilliard equations. Discrete Contin. Dyn. Syst. (2010) 28: 1669-1691.
    [17] An energy-stable and convergent finite-difference scheme for the phase field crystal equation. SIAM J. Numer. Anal. (2009) 47: 2269-2288.
    [18] Linear, first and second-order, unconditionally energy stable numerical schemes for the phase field model of homopolymer blends. J. Comput. Phys. (2016) 327: 294-316.
  • This article has been cited by:

    1. Duo Cao, Jie Shen, Jie Xu, Computing interface with quasiperiodicity, 2021, 424, 00219991, 109863, 10.1016/j.jcp.2020.109863
    2. Youngjin Hwang, Ildoo Kim, Soobin Kwak, Seokjun Ham, Sangkwon Kim, Junseok Kim, Unconditionally stable monte carlo simulation for solving the multi-dimensional Allen–Cahn equation, 2023, 31, 2688-1594, 5104, 10.3934/era.2023261
  • Reader Comments
  • © 2020 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)
通讯作者: 陈斌, bchen63@163.com
  • 1. 

    沈阳化工大学材料科学与工程学院 沈阳 110142

  1. 本站搜索
  2. 百度学术搜索
  3. 万方数据库搜索
  4. CNKI搜索

Metrics

Article views(2589) PDF downloads(179) Cited by(2)

Figures and Tables

Figures(6)  /  Tables(1)

Other Articles By Authors

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog