Cahn-Hilliard方程的一个超紧致有限差分格式.pdf

第 38 卷 第 1 期2024 年 1 月山 东 理 工 大 学 学 报(自 然 科 学 版)Journal of Shandong University of Technology(Natural Science Edition)Vol.38 No.1Jan.2024收稿日期:20221209基金项目:陕西省自然科学基金项目(2018JQ1043)第一作者:栗雪娟,女,lxj_;通信作者:王丹,女,1611182118 文章编号:1672-6197(2024)01-0073-06Cahn-Hilliard 方程的一个超紧致有限差分格式栗雪娟,王丹(西安建筑科技大学 理学院,陕西 西安 710055)摘要:研究四阶 Cahn-Hilliard 方程的数值求解方法。给出组合型超紧致差分格式,将其用于四阶 Cahn-Hilliard 方程的空间导数离散,采用四阶 Runge-Kutta 格式离散时间导数,将二者结合得到四阶 Cahn-Hilliard 方程的离散格式,并给出了该格式的误差估计。通过编程计算得到其数值解,并与精确解进行对比,结果表明本文的数值方法误差小,验证了所提方法的有效性和可行性。关键词:四阶 Cahn-Hilliard 方程;组合型超紧致差分方法;四阶 Runge-Kutta 方法;误差估计中图分类号:TB532.1;TB553文献标志码:AA supercompact finite difference scheme for Cahn-Hilliard equationsLI Xuejuan,WANG Dan(School of Science,Xian University of Architecture and Technology,Xian 710055,China)Abstract:A numerical method for solving the fourth order Cahn-Hilliard equation is studied.The combi-national ultra-compact difference scheme is given and applied to the spatial derivative discretization of thefourth order Cahn-Hilliard equation.The fourth-order Runge-Kutta scheme is used to discrete time deriv-atives.The discrete scheme of the fourth order Cahn-Hilliard equation is obtained by combining the twomethods,and the error estimate of the scheme is given.Finally,the numerical solution is obtained byprogramming and compared with the exact solution.The results show that the numerical method in thispaper has a small error,verifying the effectiveness and feasibility of the proposed method.Keywords:fourth order Cahn-Hilliard equation;combinational supercompact difference scheme;fourthorder Runge-Kutta;error estimation 本文考虑的四阶 Cahn-Hilliard 方程为ut-f u()xx+kuxxxx=0,x 0,2,t 0,u x,0()=u0 x(),x 0,2,u 0,t()=0,u 2,t()=0,t 0,(1)式中:求解区域为 0,2,且 kn 0;f u()为光滑函数;u0 x()表示 t=0 时刻的初值;ut表示 u 关于时间 t 求偏导数,ut=ut;f u()xx表示 f u()关于 x求二阶偏导数,f u()xx=2f u()x2;uxxxx表示 u 关于 x求四阶偏导数,uxxxx=4ux4;u 是混合物中某种物质的浓度,被称为相变量。1958 年,Cahn 和 Hilliard 提出 Cahn-Hilliard 方程,该方程最早被用来描述在温度降低时两种均匀的混合物所发生的相分离现象。随着学者对该方程的研究越来越深入,该方程的应用也越来越广泛,特别是在材料科学和物理学等领域中有广泛的应用1-3。Cahn-Hilliard 方程的数值解法目前已有很多研究,文献4使用了全离散有限元方法,文献5使用 了 一 类 二 阶 稳 定 的 Crank-Nicolson/Adams-Bashforth 离散化的一致性有限元逼近方法,文献6-7使用了有限元方法,文献8使用了不连续伽辽金有限元方法,文献9使用了 Cahn-Hilliard 方程的完全离散谱格式,文献10使用了高阶超紧致有限差分方法,文献11使用了高阶优化组合型紧致有限差分方法。综上所述,本文拟对 Cahn-Hilliard 方程构造一种新的超紧致差分格式,将空间组合型超紧致差分方法和修正的时间四阶 Runge-Kutta 方法相结合,求解 Cahn-Hilliard 方程的数值解,得到相对于现有广义格式精度更高的数值求解格式,并对组合型超紧致差分格式进行误差估计,最后通过数值算例验证该方法的可行性。1 高阶精度数值求解方法1.1 空间组合型超紧致差分格式早期的紧致差分格式是在 Hermite 多项式的基础上构造而来的,Hermite 多项式中连续三个节点的一阶导数、二阶导数和函数值的数值关系可以表示为1k=-1akfi+k+bkfi+k+ckfi+k()=0。(2)1998 年,Krishnan 提出如下紧致差分格式:a1fi-1+a0fi+a2fi+1+h b1fi-1+b0fi+b2fi+1()=1hc1fi-2+c2fi-1+c0fi+c3fi+1+c4fi+2(),(3)式中:h 为空间网格间距;a1,a0,a2,b1,b0,b2,c1,c2,c0,c3,c4均表示差分格式系数;fi表示 i 节点的函数值;fi和 fi分别表示 i 节点的一阶导数值和二阶导数值;fi-1,fi-2,fi+1,fi+2分别表示 i 节点依次向前两个节点和依次向后两个节点的函数值;fi-1,fi+1分别表示 i 节点依次向前一个节点和依次向后一个节点的一阶导数值;fi-1,fi+1分别表示 i 节点依次向前一个节点和依次向后一个节点的二阶导数值。式(2)对应f(x)展开以xi为邻域的泰勒级数为f x()=f xi()+hf xi()+h2f xi()2!+h3fxi()3!+h4f4()xi()4!+h5f5()xi()5!+h6f6()xi()6!+h7f7()xi()7!。(4)差分格式的各项系数由式(3)决定,可得到如下的三点六阶超紧致差分格式:716fi+1+fi-1()+fi-h16fi+1-fi-1()=1516hfi+1-fi-1(),98hfi+1-fi-1()+fi-18fi+1+fi-1()=3h2fi+1-2fi+fi-1()(5)为优化三点六阶紧致差分格式,并保持较好的数值频散,将迎风机制12引入式(5),构造出如下三点五阶迎风型超紧致差分格式:78fi-1+fi+h19fi-1-718fi-172fi+1()=1h-10148fi-1+73fi-1148fi+1(),25fi-1+fi+1h1910fi-1+165fi+910fi+1()=1h2-135fi-1-45fi+175fi+1()。(6)左右边界可达到三阶精度紧致格式:f1-132f2+f3()+3h4f2+f3()=-12hf3-f2(),f1+3728hf3-f2()+3914hf1-3356f3-f2()=f3-2f1+f2(),(7)fN-132fN-2+fN-1()-3h4fN-2+fN-1()=12hfN-2-fN-1(),fN-3728h(fN-2-fN-1)-3914hfN-3356(fN-2-fN-1)=1314h2fN-2-2fN+fN-1()。(8)上述组合型超紧致差分格式只需要相邻的三个节点便可以同时求得一阶导数和二阶导数的五阶精度近似值,比普通差分格式的节点更少,降低了计算量。为便于编程计算,将上述构造的组合型超紧致差分格式重写为矩阵表达形式。假设 U 为位移矩阵,其大小为 m n,则求一阶导数和二阶导数的离47山 东 理 工 大 学 学 报(自 然 科 学 版)2024 年散过程可以用矩阵运算表示为AF=BU,(9)结合内点的三点五阶迎风型超紧致差分格式和边界点的三点三阶差分格式,组成式(9)中等式左边的矩阵 A 和等式右边的矩阵 B,大小分别为 2m 2n和 2m n;F 为奇数行为空间一阶导数和偶数行为空间二阶导数组成的矩阵,大小为 2m n。以上矩阵分别为:A=10-13/23h/4-13/23h/439/14h1-37/28h33/5637/28h-33/567/8h/91-7h/180-h/7219/10h2/516/5h19/1007/8h/91-7h/180-h/7219/10h2/516/5h19/1007/8h/91-7h/180-h/7219/10h2/516/5h19/100-13/2-3h/4-13/2-3h/410-37/28h-33/5637/28h33/56-39/14h1,(10)F=ux()1,1ux()1,2ux()1,n-1ux()1,n2ux2()1,12ux2()1,22ux2()1,n-12ux2()1,nux()m,1ux()m,2ux()m,n-1ux()m,n2ux2()m,12ux2()m,22ux2()m,n-12ux2()m,n,(11)B=012/h-12/h-13/7h213/14h213/14h2-101/48h7/3h-11/48h-13/5h2-4/5h217/5h2-101/48h27/3h-11/48h-13/5h2-4/5h217/5h2-101/48h7/3h-11/48h-13/5h2-4/5h217/5h2012/h-12/h-13/7h213/14h213/14h2,(12)U=u1,1u1,2u1,n-1u1,nu2,1u2,2u2,n-1u2,num-1,1um-1,2um-1,n-1um-1,num,1um,2um,n-1um,n。(13)由式(9)可得F=A-1BU。(14)解线性代数方程组(9)可得 Cahn-Hilliard 方程的空间一阶导数和二阶导数。对于四阶导数,可将已求得的二阶导数替代式(14)中的 U,再次使用式(14)进行求取。57第 1 期 栗雪娟,等:Cahn-Hilliard 方程的一个超紧致有限差分格式1.2 时间离散格式在对很多偏微分方程的数值求解中不仅需要高精度的空间离散格式,同时还需要高精度的时间离散格式。普通的一阶精度时间离散格式显然满足不了高精度计算要求,因此本文选用时间四阶 Runge-Kutta 格式进行时间离散。Runge-Kutta 方法是基于欧拉方法改进后的求解偏微分方程的常用方法,这种方法不仅计算效率高,而且稳定性好。格式的推算过程如下:假设求解方程为ut+F u()=0,(15)式中 F 是对空间变量的微分算子,则修正的四阶Runge-Kutta 格式为u0i=uni,u1i=uni-t4F u()()0i,u2i=uni-t3F u()()1i,u3i=uni-t2F u()()2i,un+1i=uni-t F u()()3i。(16)1.3 误差估计以五阶精度将 fi-1,fi+1,fi-1,fi+1泰勒级数展开:fi-1=fi-hfi+h22!f(3)i-h33!f(4)i+h44!f(5)i-h55!f(6)i,fi+1=fi+hfi+h22!f(3)i+h33!f(4)i+h44!f(5)i+h55!f(6)i,fi-1=fi-hf(3)i+h22!f(4)i-h33!f(5)i+h44!f(6)i-h55!f(7)i,fi+1=fi+hf(3)i+h22!f(4)i+h33!f(5)i+h44!f(6)i+h55!f(7)i。(17)将式(17)代入式(6),所求得组合型超紧致差分格式的一阶导数及二阶导数对应的截断误差为:78fi-1+fi+h19fi-1-718fi-172fi+1()=1h-10148fi-1+73fi-1148fi+1()+78640f6()ih5,25fi-1+fi+1h1910fi-1+165fi+910fi+1()=-135fi-1-45fi+175fi+1()-5125200f7()ih5,(18)78640f6()ih5 8.101 10-4f6()ih5,5125200f7()ih5 2.023 10-3f7()ih5。(19)使用组合型超紧致差分格式的好处是在每一个网格点上存在一个一阶和二阶连续导数的多项式。本文比较了组合型超紧致差分格式和现有广义格式的一阶导数和二阶导数的截断误差:fi+fi+1+fi-1()+fi+2+fi-2()=afi+1-fi-12h+bfi+2-fi-24h+cfi+3-fi-36h,fi+fi+1+fi-1()+fi+2+fi-2()=afi+1-2fi+fi-1h2+bfi+2-2fi+fi-24h2+cfi+3-2fi+fi-39h2,(20)式中参数,a,b,c 在各种格式中取不同的值(表1,表 2)。本文发现在各种方案中,组合型超紧致差分格式的截断误差最小。表 1 不同格式一阶导数的截断误差格式abc截断误差二阶中心0010013!f3()ih2标准 Pade格式 1/403/200-15f5()ih4六阶中心003/2-3/5 1/1036 17!f7()ih6五阶迎风14316!f6()ih5表 2 不同格式二阶导数的截断误差格式abc截断误差二阶中心001002 14!f4()ih2标准 Pade格式 1/1006/50018516!f6()ih4六阶中心003/2-3/51/1072 18!f8()ih6五阶迎风16517!f7()ih567山 东 理 工 大 学 学 报(自 然 科 学 版)2024 年2 数值算例误差范数 L1和 L2的定义为:L1=1NNi=1u-U,L2=1NNi=1u-U()2。对四阶Cahn-Hilliard 取f u()=u2,k=2,在边界条件u 0,t()=u 2,t()=0 下的计算区域为 0,2,方程的精确解为 u x,t()=e-tsinx2,数值解为 U。对给出的数值算例,计算误差范数 L1和 L2,并采用四种方法进行数值模拟,对其数值结果进行误差分析和对比,结果见表 3,本文所使用方法效果最佳,由此证明所提方法的有效性和可行性。表 3 0.5 s 时刻精确度测试结果(N=10)方法L1误差L2误差间断有限元格式1.562 35 10-21.378 23 10-2普通中心差分格式1.666 67 10-18.333 33 10-2紧致差分格式7.142 86 10-31.785 71 10-3组合型超紧致差分格式6.481 48 10-36.349 21 10-4 用本文提出的式(6)式(8)和式(16)计算算例,图 1图 3 给出了不同时刻数值解与精确解的(a)精确解(b)数值解图 1 0.1 s 的精确解与数值解(a)精确解(b)数值解图 2 0.5 s 的精确解与数值解(a)精确解(b)数值解图 3 1 s 的精确解与数值解77第 1 期 栗雪娟,等:Cahn-Hilliard 方程的一个超紧致有限差分格式对比图,可以看出,数值解与精确解吻合很好,表明本文给出的数值格式是可行的,并且精度较高。3 结论本文研究了组合型超紧致差分方法和四阶Runge-Kutta 方法,并将其运用于四阶 Cahn-Hilliard方程的数值求解,通过研究与分析,得到如下结论:1)使用泰勒级数展开锁定差分格式系数,得到本文的组合型超紧致差分格式精度更高,误差更小。2)在边界点处有效地达到了降阶,并提高了精度。3)通过数值算例验证了数值格式的有效性。4)预估该方法可应用于高阶偏微分方程的数值求解。参考文献:1HUANG Q M,YANG J X.Linear and energy-stable method with en-hanced consistency for the incompressible Cahn-Hilliard-Navier-Stokestwo-phaseflowmodelJ.Mathematics,2022,10(24):4711.2AKRIVIS G,LI B Y,LI D F.Energy-decaying extrapolated RK-SAVmethods for the allen-Cahn and Cahn-Hilliard equationsJ.SIAMJournal on Scientific Computing,2019,41(6):3703-3727.3YOUNAS U,REZAZADEH H,REN J,et al.Propagation of diverseexact solitary wave solutions in separation phase of iron(Fe-Cr-X(X=Mo,Cu)for the ternary alloysJ.International Journal of Mod-ern Physics B,2022,36(4):2250039.4 HE R J,CHEN Z X,FENG X L.Error estimates of fully discretefinite element solution for the 2D Cahn-Hilliard equation with infinitetime horizon J.Numerical Methods for Partial Differential Equati-ions,2017,33(3):742-762.5HE Y N,FENG X L.Uniform H2-regularity of solution for the 2DNavier-Stokes/Cahn-Hilliard phase field model J.Journal of Math-ematical Analysis and Applications,2016,441(2):815-829.6WEN J,HE Y N,HE Y L.Semi-implicit,unconditionally energy sta-ble,stabilized finite element method based on multiscale enrichmentfortheCahn-Hilliard-Navier-Stokesphase-fieldmodel J.Computers and Mathematics with Applications,2022,126:172-181.7 MESFORUSH A,LARSSON S.A posteriori error analysis for theCahn-Hilliard equation J.Journal of Mathematical Modeling,2022,10(4):437-452.8XIA Y,XU Y,SHU C W.Local discontinuous Galerkin methods forthe Cahn-Hilliard type equation J.Journal of Computational Phys-ics,2007,227(1):472-491.9CHEN L,L S J.A fully discrete spectral scheme for time fractionalCahn-Hilliard equation with initial singularityJ.Computers andMathematics with Applications,2022,127:213-224.10周诚尧,汪勇,桂志先,等.二维黏弹介质五点八阶超紧致有限差分声波方程数值模拟J.科学技术与工程,2020,20(1):54-63.11汪勇,徐佑德,高刚,等.二维黏滞声波方程的优化组合型紧致有限差分数值模拟J.石油地球物理勘探,2018,53(6):1152-1164,1110.12程晓晗,封建湖,郑素佩.求解对流扩散方程的低耗散中心迎风格式J.应用数学,2017,30(2):344-349.(编辑:杜清玲)87山 东 理 工 大 学 学 报(自 然 科 学 版)2024 年