基于局部加密纯无网格法非线性Cahn-Hilliard方程的模拟*

任金莲 蒋戎戎 陆伟刚 蒋涛

(扬州大学数学科学学院,扬州 225002)

(扬州大学水利科学与工程学院,扬州 225002)

为数值求解描述不同物质间相位分离现象的高阶非线性Cahn-Hilliard(C-H)方程,发展了一种基于局部加密纯无网格有限点集法(local refinement finite pointset method,LR-FPM).其构造过程为: 1)将 C-H 方程中四阶导数降阶为两个二阶导数,连续应用基于Taylor展开和加权最小二乘法的FPM离散空间导数;2)对区域进行局部加密和采用五次样条核函数以提高数值精度;3)局部线性方程组求解中准确施加含高阶导数Neumann边值条件.随后,运用LR-FPM求解有解析解的一维/二维 C-H方程,分析粒子均匀分布/非均匀分布以及局部粒子加密情况的误差和收敛阶,展示了LR-FPM较网格类算法在非均匀布点情况下的优点.最后,采用LR-FPM对无解析解的一维/二维 C-H方程进行了数值预测,并与有限差分结果相比较.数值结果表明,LR-FPM方法具有较高的数值精度和收敛阶,比有限差分法更易数值实现,能够准确展现不同类型材料间相位分离非线性扩散现象随时间的演化过程.

1 引 言

热力学、流体力学和生物数学等领域中常涉及复杂相位分离现象的非线性扩散问题[1−4],这类问题常用含四阶导数 Cahn-Hilliard(C-H)方程[3−6]来描述.然而,许多情况下C-H方程含有的非线性项使其难以得到其光滑理论解[6−10].C-H方程的数值算法及模拟研究已成为科学计算领域中一个国际热点问题.目前,已提出多种数值方法对C-H方程或对高维情况下非线性扩散过程进行求解或数值预测,如有限差分法、有限元法、蒙特卡罗法、谱方法、修正欧拉法和无单元伽辽金法等[6−22].上述方法均基于网格或背景网格,在非均匀分布节点或多维复杂区域及含高阶导数边界条件情况下程序实现比较复杂,算法的灵活推广应用性上受到很大限制.因此,近年来纯无网格粒子方法,如光滑粒子动力学方法(smoothed particle hydrodynamics,SPH)和有限点集法(finite pointset method,FPM)等[23−31],以其不受网格限制和易推广应用到复杂高维区域问题上的优势受到普遍关注,亦在偏微分方程的数值求解中得到了许多应用.值得注意的是,上述纯无网格粒子法对含高阶导数多维偏微分方程的研究还处在探索阶段,特别是国际上鲜有关于C-H方程的纯无网格粒子法研究的报道.

纯无网格FPM法在含高阶导数C-H方程的数值模拟上还处在试探性研究阶段,这与该方法的数值精度和稳定性方面有待进一步研究有关,其数值精度的提高和其对非线性C-H方程的数值模拟均是两个研究热点.FPM法数值模拟C-H方程较网格类方法的优点主要在于: 易处理非均匀节点分布或局部加密的情况、易离散求解高阶空间导数且不依赖于网格、易准确施加含高阶导数Neumann边界.

基于上述原因,本文结合高阶C-H方程的特点,提出一种能够有效、准确地模拟非线性C-H方程局部加密FPM方法(local refinement finite pointset method,LR-FPM): 首先将四阶空间导数分解为两个二阶导数,区域离散时采用局部加密;其次采用五次样条核函数,基于Taylor展开和加权最小二乘思想的FPM法对两个二阶空间导数依次离散,并将含高阶导数的Neumann边界条件准确施加到离散格式中.通过数值算例展现了提出的LR-FPM方法能准确地求解一维/二维非线性C-H方程,具有二阶收敛性,能够准确可靠地预测C-H方程描述的非线性热扩散现象随时间演化过程.

2 C-H方程

C-H方程[3−6,11]是一个含四阶导数的非线性偏微分方程,常用来描述一种旋量分解的相位分离现象或工程中的非线性扩散现象.C-H方程包含非线性项、四阶导数和Neumann边值条件中三阶导数,使其成为验证一种数值模拟算法能力的挑战性算例[6,9−11].本文考虑如下形式的C-H方程:

其中 u(x,t)是二元混合物中组分的浓度,F′(u(x,t))是自由能函数,ε为梯度界面能系数,Mr为恒定的迁移率.根据文献[11],方程(1)常被分裂为含二阶空间导数的两个偏微分方程:

其中 μ 是局部化学势,对应的初边值条件为:

初始条件

齐次Neumann边界条件

或周期边界条件

其中∇ u和∇μ为梯度向量,n 为 ∂ Ω上的外法向量.

C-H方程对应Ginzburg–Landau自由能方程为:

3 局部加密FPM离散方法(LR-FPM)

本文主要考虑含高阶导数非线性问题的粒子方法模拟,基于高阶导数降阶和局部加密思想,拓展FPM法得到一种能够准确模拟C-H方程的LR-FPM.给出的LR-FPM离散算法的基本思想概括为: a)基于文献[11]将含四阶空间导数(1)式分裂为两个含二阶空间导数的(2)式,对(2)式两个方程中空间导数依次拓展应用FPM离散格式;b)FPM离散中采用较传统高斯核函数具有较高精度的五次样条核函数[23],两个Neumann边值条件(4)式也依次准确地施加到形成的线性方程组中;c)时间导数上采用二阶精度预估校正格式,对计算区域进行局部加密以提高数值精度的同时还降低了计算量.

3.1 FPM离散格式

对四阶导数降阶后得到的两个二阶导数方程(2)的空间导数离散求解,连续两次应用FPM,函数一阶/二阶导数的显式FPM离散近似思想基于Taylor展开和加权最小二乘法(详见文献[28,29,31]).

设 u(x)为空间函数,计算区域 Ω ⊂Rd(d=1,2,3)内可以不依赖网格任意布点 xj,j=1,2,···,N,(N为Ω中的粒子总数),函数在 x 处的一阶/二阶导数以加权函数支持域内n个相邻粒子xi(i=1,2,···,n)来近似.加权最小二乘近似中采用较传统高斯核函数[28−30]具有较高精度的五次样条核函数[23],其形式如下:

其中 q=rij/h,rij=|xi-xj|,αd是正常数,αd在一维、二维空间中分别为 1 20/h,7/(478πh2),h 为光滑长度,此处取 h ≈1.1d0(d0为分布粒子的初始距离).在粒子均匀分布情况下的支持域范围是以3h为半径的圆,非均匀粒子分布或局部加密时则采用平均方式的变光滑长度(详见文献[23]).

考虑相邻粒子 xi在x处的Taylor 展开,可得

其中 ei是Taylor展开式的误差余项,符号x1i,x2i分别用来表示点 xi的 x,y 分量,k,l 表示分量的偏导数.(9)式可化为

其中

dxi,dyi分别表示xi-x,yi-y,(i=1,2,···,n).

对于未知函数 u 的一阶/二阶导数通过误差ei加权最小二乘法来计算,可得

(11)式可以写成J=(Ma-b)TW(Ma-b),其中 W 为对角矩阵,对角线元素为 w1,w2,· ··,wn.

根据 J 的极小值原理,得到

(12)式涉及 5 ×5 局部系数矩阵,可通过其求解得到一阶/二阶导数近似值.若求解带有Neumann边值条件的问题,上述涉及的矩阵 M,W和向量b,e 均要发生变化(详见 3.2 节).

3.2 边界处理及粒子分布局部加密

众所周知,偏微分方程数值计算中边值条件的处理对数值模拟的精度和稳定性至关重要,也是计算方法对其准确处理的难点之一.本文对方程(2)的数值模拟中涉及三种初边值条件的处理,对于初始条件(3)式和周期边界条件(5)式可以根据文献[7,23]准确的施加.对Neumann边界条件,离散并准确施加到方程(12)中的过程是:a)在(2)式的第一次FPM离散过程中将Neumann条件中∇u·n离散为0=uxnx+uyny(nx,ny是单位外法向量),将其与(9)式进行联合,此时式中的矩阵M增加一行元素为(nx,ny,0,0,0),W对角线上元素变为w1,w2,···,wn,1,b=(u1-u,u2-u,u3-u,···,un-u,0)T,e=(e1,e2,e3,···,en,en+1)T;b)在(2)式的第二次FPM离散过程中,Neumann条件中的离 散为然后将其与(9)式联合,求解局部线性方程组时矩阵 M,W和向量 b,e 的变化与a)类似.第4节数值验证了上述边值条件的处理是准确的,特别对于Neumann边界条件,且较网格类有限差分法容易施加.

完全不依赖网格的纯无网格FPM粒子法,在离散格式构造中可在计算区域内任意布置粒子,使其在粒子分布非均匀情况下的编程计算实现比网格类方法容易,且易推广应用到复杂高维非规则区域问题的模拟.为体现给出的LR-FPM模拟C-H方程的上述优点,数值模拟中考虑了粒子非均匀分布情况,为了兼顾提高计算精度和降低计算量对计算区域进行局部加密.值得注意的是粒子分布局部加密情况下,在粗粒子分布和加密粒子分布相交区域里采用加权平均的变光滑长度(详见文献[23]).由于边界处支持域内粒子缺失会降低数值精度和稳定性,以及模拟区域中高峰值附近数值近似对整个计算精度有重要影响,因此,数值算例中主要针对边界处或高峰值函数附近小区域进行粒子局部加密(见图3(b)).

值得注意的是,在第4节和第5节的数值算例模拟中,均采用FPM离散(12)式结合3.2节边界处理和时间项二阶精度离散格式对C-H方程(1)—(4)式进行离散求解.

4 FPM方法精度和收敛性分析

本节数值求解有解析解但边界条件不同的一/二维C-H方程,并结合数值结果分析LR-FPM方法的数值误差和收敛速度.所模拟的问题是一维带周期边值、二维带Neumann边值问题,并与解析解做比较.为体现本文给出的方法相对于网格类方法的优点,着重讨论了在FPM方法的基础上进行局部加密分布和非均匀分布情况.局部加密采用了疏密点相结合的分布方式,这种处理方式的优点在于: 1)相对于全局加密,局部加密增加数值精度的同时减少了计算量;2)相对于网格类方法,局部加密简单易行.

为了分析LR-FPM方法的精度和收敛性,定义如下误差(精确解与数值解L2-范数误差)和收敛阶:

其中 u 为数值解;U 为解析解;d1,d2为不同的粒子初始间距.

4.1 一维周期边界C-H方程

考虑区域 Ω =[0,2π] 上带周期边值条件的C-H方程[7],其对应的方程和初边值条件分别为

该算例模拟中,均匀分布情况下粒子初始间距为 d1=π/32,时间步长为 d t=10-4.局部加密情况下,在 [ 0 ,10π/32],[ 5 4π/32,2π] 处加密,加密区域的粒子间距为初始间距的1/2,即加密区域间距为d2=π/64,[ 1 0π/32,54π/32] 内仍采用粒子初始间距d1,时间步长为 d t=10-4,加密临界点光滑长度取h≈0.5×(1.1×d1+1.1×d2).图 1展示了不同时刻下均匀分布和局部加密两种情况下的FPM数值解,并与解析解作对比.可知两种情况下的LRFPM结果均与解析解吻合.

图1 几个不同时刻下均匀分布、局部加密情况下的数值解和解析解Fig.1.Comparisons between the present numerical results and analytical solutions with different times under the uniform and local refinement particle distributions.

为体现提出的FPM方法模拟周期性问题的数值精度和收敛速度,图2给出了较短时间内不同粒子数下数值误差-时间曲线,L2-范数误差 E2随着粒子数的增加而减小,粒子数较多情况下具有更好的数值收敛性.值得注意的是在粒子数N=129时L2-范数误差 E2先减小后增大,这是由空间步长较小和随时间延长函数值变小情况下计算机的舍入误差导致,但随时间延长中误差 E2基本稳定在同一量级,仍符合误差随粒子增加而减小的数值计算理论.同时为了进一步体现LR-FPM模拟周期性问题的数值精度,表1列出了不同粒子间距下t=0.5s时的L2-范数误差和收敛阶.由表1可知,LR-FPM具有二阶收敛速度,与文献[7]中有限差分方法的收敛阶基本一致,从而表明LR-FPM模拟C-H方程时具有较高的数值精度,较网格类方法具有更好的灵活性和推广应用性.

图2 不同粒子数下的收敛速度随时间的变化Fig.2.The numerical convergence versus time under different particle numbers.

表1 t=0.5s 时不同粒子间距情况下的 L2-范数误差 E2和收敛阶Table 1.The L2-norm error E2 and convergence rate at t=0.5s .

表2列出了不同时刻均匀分布与局部加密情况的L2-范数误差 E2,由表2可知局部加密情况下的误差远远小于均匀分布情况下的误差.这表明LR-FPM不仅具有更好的灵活推广应用性,而且能保持较高的精度.

表2 不同时刻下粒子均匀分布与局部加密情况下的L2-范数误差 E2 对比Table 2.The L2-norm error E2 at different times under the uniform and local refinement particle distributions.

4.2 二维Neumann边界C-H方程

为体现提出的LR-FPM方法求解带第二类边界非线性C-H方程的准确性,本节考虑正方形区域 Ω =[0,1]×[0,1] 的第二类边界条件的非线性C-H方程,其对应的方程和初边值条件[6]为

图3展示了本算例模拟中采用的三种粒子分布方式,其中均匀分布方式是分别沿 x,y 方向每隔0.04的距离分布一个粒子,局部加密分布采用四角和中间加密,分别在 [0,0.1] × [0,0.1],[0,0.1] ×[0.9,1.0],[0.9,1.0] × [0,0.1],[0.9,1.0] × [0.9,1.0],以及 [0.4,0.6] × [0.4,0.6]采用加密形式,而非均匀分布以(0.5,0.5)为圆心,向外圆形分布粒子,相邻的两个圆距离相等,其余区域上仍采用均匀布点的方式.

图3 不同粒子分布(a)粒子均匀分布;(b)粒子局部加密分布;(c)粒子非均匀分布Fig.3.Different cases of particle distributions:(a)Uniform case;(b)local refinement case;(c)non-uniform case.

该算例模拟中,均匀分布采用粒子初始间距为d1=0.04,时间步长为 d t=10-6,局部加密采用加密区域粒子间距为 d2=0.025,其余区域粒子间距仍为 d1=0.04,加密临界点光滑长度取h≈0.5×(1.1×d1+1.1×d2).

图4展示了 t=0.01s 时刻均匀分布和局部加密两种粒子分布情况下 x=y 处的数值解,并与解析解进行对比.可以看出,随时间演化,两种情况下的数值结果都与解析解吻合,但局部加密情况下的结果更接近解析解.

图4 均匀分布与局部加密情况下的数值结果Fig.4.The present numerical results under the uniform and local refinement particle distributions.

表3列出了粒子间距为 d0=0.04 在FPM方法下采用五次样条核函数与高斯核函数情况下的L2-范数误差 E2对比.由表 3 可知,采用五次样条核函数的误差比采用高斯核函数的误差小.因此,为提高数值精度,模拟中均采用五次样条核函数.

表3 初始间距 d0=0.04 情况下五次样条核函数与高斯核函数的L2-范数误差 E2 对比Table 3.The L2-norm error with quintic spline kernel and gaussian kernel functions at d0=0.04 .

为进一步体现提出的LR-FPM方法模拟Neumann边界问题的数值精度和收敛速度,表4列出了模拟较短时间内数值结果的误差和收敛阶,通过表4可得: 数值误差随着粒子数的增加而减小,且 d0=1/20 到 d0=1/40的误差变化比d0=1/40 到 d0=1/60 的误差变化大.给出的LRFPM方法模拟Neumann边界条件下C-H方程具有二阶收敛速度.

表5列出了粒子均匀分布、局部加密分布与非均匀分布情况下的L2-范数误差 E2.通过观察表5,本文给出的粒子方法在粒子分布均匀和非均匀情况下得到的数值误差接近,局部加密情况下的误差远远小于均匀分布情况下的误差,表明局部加密能有效地提高数值精度.

表4 t=0.01s 时刻下不同粒子间距的 L2-范数误差 E2和收敛阶Table 4.The L2-norm error E2 and convergence rate at t=0.01s .

表5 粒子均匀分布、局部加密分布与非均匀分布情况下的L2-范数误差 E2 对比Table 5.The L2-norm error E2 at different times under the uniform,local refinement,and non-uniform particle distributions.

为进一步体现粒子方法在粒子非均匀分布情况下数值模拟精度,表6列出了 t=0.01s 时不同粒子间距非均匀分布情况下的误差和收敛阶.由表4、表5和表6可知,粒子方法在均匀分布和非均匀分布情况下得到的数值结果误差都比较接近,表明了该方法易推广应用到非均匀区域问题的模拟,且保持较好的精度.本文的粒子方法不受网格限制,可以在模拟区域任意粒子布置情况下对问题进行模拟实现,较网格类方法更容易被推广应用于复杂非规则区域上C-H问题的模拟.

表6 t=0.01 s 时不同粒子间距非均匀分布情况下的L2-范数误差 E2和收敛阶Table 6.The L2-norm error E2 and convergence rate at t=0.01 s under non-uniform particle distribution.

5 C-H 方程数值模拟

为了更好地体现提出的LR-FPM粒子方法数值模拟C-H问题的准确性,本节分别选取一维/二维Neumann边界无解析解C-H问题,并与文献[8]中基于离散变分导数的有限差分(finite difference method,FDM)法模拟结果做对比.

5.1 一维Neumann边界C-H方程

考虑区域 Ω=[-1,1] 上带有Neumann边值条件的C-H方程,其表达式为

对应的初值条件为u0(x)=0.53x+0.47sin(-1.5πx)和边值条件为

该算例为一维无解析解的带Neumann边界C-H方程,它将描述复杂的相位分离现象,常被用来证明其保持方程的质量守恒性质和能量耗散性质.本节运用LR-FPM方法对 ε2=0.3 情况下算例进行了模拟,并与文献[8]中FDM方法结果进行对比.图 5 给出了 ε2=0.3,时间步长为 d t=10-6,初始时刻与t=0.2s,t=2s 三个不同时刻两种数值方法的模拟结果.由图5可知,LR-FPM方法得到的带Neumann边界C-H方程的数值解与FDM 结果吻合.图6给出了ε2=0.03,t=0.2s 时刻下粒子局部加密情况,在 [–1.0,–0.8],[0.8,1.0]处加密,每隔 0.02 的距离布置粒子,在(–0.8,0.8)内每隔0.04的距离布置粒子.由图6可知,粒子均匀分布与局部加密情况下的数值结果均与FDM结果吻合,但局部加密情况下的模拟结果更接近于FDM的结果.这进一步表明提出的LR-FPM粒子方法模拟C-H方程是可靠的.

图5 ε2=0.3 时不同时刻 FDM 结果与 LR-FPM 结果Fig.5.The numerical results obtained using FDM and LRFPM at different times with ε2=0.3 .

图6 ε2=0.03,t=0.2s 时刻下均匀分布与局部加密情况下数值结果对比Fig.6.The present numerical results under uniform and local refinement particle distributions at ε2=0.03,t=0.2 s.

5.2 二维Neumann边界C-H方程

为体现提出的FPM方法求解带Neumann边界非线性C-H方程的准确性,本节考虑正方形区域 Ω =[0,0.8]×[0,0.8] 的第二类边界条件的非线性C-H方程,其对应的方程和初边值条件[11]为

初值条件为

和Neumann边界条件

图7 t=0.1s 时刻 FPM 方法模拟结果Fig.7.The FPM result at t=0.1s .

图7给出了 t=0.1s 时刻LR-FPM方法的模拟结果,与文献[11]中数值结果相比较,可以看出两种方法的结果接近.图 8 给出了 t=0.08 s时刻GRBF方法[11]等值线图(图8(a))和本文方法在三种粒子分布情况下的数值等值线分布图(图8(b)—图8(d)).通过观察图8可知,本文方法在粒子均匀分布、局部加密以及非均匀分布情况下得到的带Neumann边界C-H方程的数值等值线分布均与GRBF结果吻合,表明本文方法数值预测带高阶导数Neumann边界C-H问题是可靠的.

图9给出了 t=0.08s 时刻下不同粒子间距均匀分布的情况以及局部加密分布下的数值收敛性.由图9可知,LR-FPM数值预测二维C-H问题时是收敛的.因此,LR-FPM能够有效、可靠地模拟C-H方程,容易实施局部加密或非均匀分布情况下的数值模拟.

图9 t=0.08s 时刻粒子局部加密分布情况下的数值收敛Fig.9.The numerical convergence obtained using the present method under different particle distributions at t=0.08s.

6 结 论

本文基于高阶导数分解和区域离散局部加密思想,首次拓展应用FPM方法对描述相位分离现象的含四阶导数C-H方程进行数值模拟研究,为复杂区域上含高阶导数非线性问题的数值预测提供了一种准确、有效的LR-FPM方法.为了验证LR-FPM方法的优势,首先采用LR-FPM求解了均匀、局部加密和非均匀离散情况下带有解析解的一维/二维C-H方程,并分析了LR-FPM的数值误差和收敛性;其次,运用LR-FPM法对无解析解C-H方程进行了数值预测,并与有限差分法结果作对比.所有数值结果表明:

1)本文给出的带有高阶导数边值条件的施加处理是准确的,且给出的LR-FPM粒子方法模拟一维/二维C-H方程具有较好二阶收敛速率;

2)基于局部加密下LR-FPM能有效地降低数值模拟误差,且在粒子分布非均匀情况下比有限差分法有更好灵活应用性;

3)给出的LR-FPM法能准确预测不同物质间相位分离非线性现象随时间的演化过程.