任金蓮 蔣戎戎 陸偉剛 蔣濤

(揚(yáng)州大學(xué)數(shù)學(xué)科學(xué)學(xué)院,揚(yáng)州 225002)

(揚(yáng)州大學(xué)水利科學(xué)與工程學(xué)院,揚(yáng)州 225002)

為數(shù)值求解描述不同物質(zhì)間相位分離現(xiàn)象的高階非線性Cahn-Hilliard(C-H)方程,發(fā)展了一種基于局部加密純無(wú)網(wǎng)格有限點(diǎn)集法(local refinement finite pointset method,LR-FPM).其構(gòu)造過(guò)程為: 1)將 C-H 方程中四階導(dǎo)數(shù)降階為兩個(gè)二階導(dǎo)數(shù),連續(xù)應(yīng)用基于Taylor展開(kāi)和加權(quán)最小二乘法的FPM離散空間導(dǎo)數(shù);2)對(duì)區(qū)域進(jìn)行局部加密和采用五次樣條核函數(shù)以提高數(shù)值精度;3)局部線性方程組求解中準(zhǔn)確施加含高階導(dǎo)數(shù)Neumann邊值條件.隨后,運(yùn)用LR-FPM求解有解析解的一維/二維 C-H方程,分析粒子均勻分布/非均勻分布以及局部粒子加密情況的誤差和收斂階,展示了LR-FPM較網(wǎng)格類(lèi)算法在非均勻布點(diǎn)情況下的優(yōu)點(diǎn).最后,采用LR-FPM對(duì)無(wú)解析解的一維/二維 C-H方程進(jìn)行了數(shù)值預(yù)測(cè),并與有限差分結(jié)果相比較.數(shù)值結(jié)果表明,LR-FPM方法具有較高的數(shù)值精度和收斂階,比有限差分法更易數(shù)值實(shí)現(xiàn),能夠準(zhǔn)確展現(xiàn)不同類(lèi)型材料間相位分離非線性擴(kuò)散現(xiàn)象隨時(shí)間的演化過(guò)程.

1 引 言

熱力學(xué)、流體力學(xué)和生物數(shù)學(xué)等領(lǐng)域中常涉及復(fù)雜相位分離現(xiàn)象的非線性擴(kuò)散問(wèn)題[1?4],這類(lèi)問(wèn)題常用含四階導(dǎo)數(shù) Cahn-Hilliard(C-H)方程[3?6]來(lái)描述.然而,許多情況下C-H方程含有的非線性項(xiàng)使其難以得到其光滑理論解[6?10].C-H方程的數(shù)值算法及模擬研究已成為科學(xué)計(jì)算領(lǐng)域中一個(gè)國(guó)際熱點(diǎn)問(wèn)題.目前,已提出多種數(shù)值方法對(duì)C-H方程或?qū)Ω呔S情況下非線性擴(kuò)散過(guò)程進(jìn)行求解或數(shù)值預(yù)測(cè),如有限差分法、有限元法、蒙特卡羅法、譜方法、修正歐拉法和無(wú)單元伽遼金法等[6?22].上述方法均基于網(wǎng)格或背景網(wǎng)格,在非均勻分布節(jié)點(diǎn)或多維復(fù)雜區(qū)域及含高階導(dǎo)數(shù)邊界條件情況下程序?qū)崿F(xiàn)比較復(fù)雜,算法的靈活推廣應(yīng)用性上受到很大限制.因此,近年來(lái)純無(wú)網(wǎng)格粒子方法,如光滑粒子動(dòng)力學(xué)方法(smoothed particle hydrodynamics,SPH)和有限點(diǎn)集法(finite pointset method,FPM)等[23?31],以其不受網(wǎng)格限制和易推廣應(yīng)用到復(fù)雜高維區(qū)域問(wèn)題上的優(yōu)勢(shì)受到普遍關(guān)注,亦在偏微分方程的數(shù)值求解中得到了許多應(yīng)用.值得注意的是,上述純無(wú)網(wǎng)格粒子法對(duì)含高階導(dǎo)數(shù)多維偏微分方程的研究還處在探索階段,特別是國(guó)際上鮮有關(guān)于C-H方程的純無(wú)網(wǎng)格粒子法研究的報(bào)道.

純無(wú)網(wǎng)格FPM法在含高階導(dǎo)數(shù)C-H方程的數(shù)值模擬上還處在試探性研究階段,這與該方法的數(shù)值精度和穩(wěn)定性方面有待進(jìn)一步研究有關(guān),其數(shù)值精度的提高和其對(duì)非線性C-H方程的數(shù)值模擬均是兩個(gè)研究熱點(diǎn).FPM法數(shù)值模擬C-H方程較網(wǎng)格類(lèi)方法的優(yōu)點(diǎn)主要在于: 易處理非均勻節(jié)點(diǎn)分布或局部加密的情況、易離散求解高階空間導(dǎo)數(shù)且不依賴(lài)于網(wǎng)格、易準(zhǔn)確施加含高階導(dǎo)數(shù)Neumann邊界.

基于上述原因,本文結(jié)合高階C-H方程的特點(diǎn),提出一種能夠有效、準(zhǔn)確地模擬非線性C-H方程局部加密FPM方法(local refinement finite pointset method,LR-FPM): 首先將四階空間導(dǎo)數(shù)分解為兩個(gè)二階導(dǎo)數(shù),區(qū)域離散時(shí)采用局部加密;其次采用五次樣條核函數(shù),基于Taylor展開(kāi)和加權(quán)最小二乘思想的FPM法對(duì)兩個(gè)二階空間導(dǎo)數(shù)依次離散,并將含高階導(dǎo)數(shù)的Neumann邊界條件準(zhǔn)確施加到離散格式中.通過(guò)數(shù)值算例展現(xiàn)了提出的LR-FPM方法能準(zhǔn)確地求解一維/二維非線性C-H方程,具有二階收斂性,能夠準(zhǔn)確可靠地預(yù)測(cè)C-H方程描述的非線性熱擴(kuò)散現(xiàn)象隨時(shí)間演化過(guò)程.

2 C-H方程

C-H方程[3?6,11]是一個(gè)含四階導(dǎo)數(shù)的非線性偏微分方程,常用來(lái)描述一種旋量分解的相位分離現(xiàn)象或工程中的非線性擴(kuò)散現(xiàn)象.C-H方程包含非線性項(xiàng)、四階導(dǎo)數(shù)和Neumann邊值條件中三階導(dǎo)數(shù),使其成為驗(yàn)證一種數(shù)值模擬算法能力的挑戰(zhàn)性算例[6,9?11].本文考慮如下形式的C-H方程:

其中 u(x,t)是二元混合物中組分的濃度,F′(u(x,t))是自由能函數(shù),ε為梯度界面能系數(shù),Mr為恒定的遷移率.根據(jù)文獻(xiàn)[11],方程(1)常被分裂為含二階空間導(dǎo)數(shù)的兩個(gè)偏微分方程:

其中 μ 是局部化學(xué)勢(shì),對(duì)應(yīng)的初邊值條件為:

初始條件

齊次Neumann邊界條件

或周期邊界條件

其中? u和?μ為梯度向量,n 為 ? Ω上的外法向量.

C-H方程對(duì)應(yīng)Ginzburg–Landau自由能方程為:

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

本文主要考慮含高階導(dǎo)數(shù)非線性問(wèn)題的粒子方法模擬,基于高階導(dǎo)數(shù)降階和局部加密思想,拓展FPM法得到一種能夠準(zhǔn)確模擬C-H方程的LR-FPM.給出的LR-FPM離散算法的基本思想概括為: a)基于文獻(xiàn)[11]將含四階空間導(dǎo)數(shù)(1)式分裂為兩個(gè)含二階空間導(dǎo)數(shù)的(2)式,對(duì)(2)式兩個(gè)方程中空間導(dǎo)數(shù)依次拓展應(yīng)用FPM離散格式;b)FPM離散中采用較傳統(tǒng)高斯核函數(shù)具有較高精度的五次樣條核函數(shù)[23],兩個(gè)Neumann邊值條件(4)式也依次準(zhǔn)確地施加到形成的線性方程組中;c)時(shí)間導(dǎo)數(shù)上采用二階精度預(yù)估校正格式,對(duì)計(jì)算區(qū)域進(jìn)行局部加密以提高數(shù)值精度的同時(shí)還降低了計(jì)算量.

3.1 FPM離散格式

對(duì)四階導(dǎo)數(shù)降階后得到的兩個(gè)二階導(dǎo)數(shù)方程(2)的空間導(dǎo)數(shù)離散求解,連續(xù)兩次應(yīng)用FPM,函數(shù)一階/二階導(dǎo)數(shù)的顯式FPM離散近似思想基于Taylor展開(kāi)和加權(quán)最小二乘法(詳見(jiàn)文獻(xiàn)[28,29,31]).

設(shè) u(x)為空間函數(shù),計(jì)算區(qū)域 Ω ?Rd(d=1,2,3)內(nèi)可以不依賴(lài)網(wǎng)格任意布點(diǎn) xj,j=1,2,···,N,(N為Ω中的粒子總數(shù)),函數(shù)在 x 處的一階/二階導(dǎo)數(shù)以加權(quán)函數(shù)支持域內(nèi)n個(gè)相鄰粒子xi(i=1,2,···,n)來(lái)近似.加權(quán)最小二乘近似中采用較傳統(tǒng)高斯核函數(shù)[28?30]具有較高精度的五次樣條核函數(shù)[23],其形式如下:

其中 q=rij/h,rij=|xi-xj|,αd是正常數(shù),αd在一維、二維空間中分別為 1 20/h,7/(478πh2),h 為光滑長(zhǎng)度,此處取 h ≈1.1d0(d0為分布粒子的初始距離).在粒子均勻分布情況下的支持域范圍是以3h為半徑的圓,非均勻粒子分布或局部加密時(shí)則采用平均方式的變光滑長(zhǎng)度(詳見(jiàn)文獻(xiàn)[23]).

考慮相鄰粒子 xi在x處的Taylor 展開(kāi),可得

其中 ei是Taylor展開(kāi)式的誤差余項(xiàng),符號(hào)x1i,x2i分別用來(lái)表示點(diǎn) xi的 x,y 分量,k,l 表示分量的偏導(dǎo)數(shù).(9)式可化為

其中

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

對(duì)于未知函數(shù) u 的一階/二階導(dǎo)數(shù)通過(guò)誤差ei加權(quán)最小二乘法來(lái)計(jì)算,可得

(11)式可以寫(xiě)成J=(Ma-b)TW(Ma-b),其中 W 為對(duì)角矩陣,對(duì)角線元素為 w1,w2,· ··,wn.

根據(jù) J 的極小值原理,得到

(12)式涉及 5 ×5 局部系數(shù)矩陣,可通過(guò)其求解得到一階/二階導(dǎo)數(shù)近似值.若求解帶有Neumann邊值條件的問(wèn)題,上述涉及的矩陣 M,W和向量b,e 均要發(fā)生變化(詳見(jiàn) 3.2 節(jié)).

3.2 邊界處理及粒子分布局部加密

眾所周知,偏微分方程數(shù)值計(jì)算中邊值條件的處理對(duì)數(shù)值模擬的精度和穩(wěn)定性至關(guān)重要,也是計(jì)算方法對(duì)其準(zhǔn)確處理的難點(diǎn)之一.本文對(duì)方程(2)的數(shù)值模擬中涉及三種初邊值條件的處理,對(duì)于初始條件(3)式和周期邊界條件(5)式可以根據(jù)文獻(xiàn)[7,23]準(zhǔn)確的施加.對(duì)Neumann邊界條件,離散并準(zhǔn)確施加到方程(12)中的過(guò)程是:a)在(2)式的第一次FPM離散過(guò)程中將Neumann條件中?u·n離散為0=uxnx+uyny(nx,ny是單位外法向量),將其與(9)式進(jìn)行聯(lián)合,此時(shí)式中的矩陣M增加一行元素為(nx,ny,0,0,0),W對(duì)角線上元素變?yōu)閣1,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離散過(guò)程中,Neumann條件中的離 散為然后將其與(9)式聯(lián)合,求解局部線性方程組時(shí)矩陣 M,W和向量 b,e 的變化與a)類(lèi)似.第4節(jié)數(shù)值驗(yàn)證了上述邊值條件的處理是準(zhǔn)確的,特別對(duì)于Neumann邊界條件,且較網(wǎng)格類(lèi)有限差分法容易施加.

完全不依賴(lài)網(wǎng)格的純無(wú)網(wǎng)格FPM粒子法,在離散格式構(gòu)造中可在計(jì)算區(qū)域內(nèi)任意布置粒子,使其在粒子分布非均勻情況下的編程計(jì)算實(shí)現(xiàn)比網(wǎng)格類(lèi)方法容易,且易推廣應(yīng)用到復(fù)雜高維非規(guī)則區(qū)域問(wèn)題的模擬.為體現(xiàn)給出的LR-FPM模擬C-H方程的上述優(yōu)點(diǎn),數(shù)值模擬中考慮了粒子非均勻分布情況,為了兼顧提高計(jì)算精度和降低計(jì)算量對(duì)計(jì)算區(qū)域進(jìn)行局部加密.值得注意的是粒子分布局部加密情況下,在粗粒子分布和加密粒子分布相交區(qū)域里采用加權(quán)平均的變光滑長(zhǎng)度(詳見(jiàn)文獻(xiàn)[23]).由于邊界處支持域內(nèi)粒子缺失會(huì)降低數(shù)值精度和穩(wěn)定性,以及模擬區(qū)域中高峰值附近數(shù)值近似對(duì)整個(gè)計(jì)算精度有重要影響,因此,數(shù)值算例中主要針對(duì)邊界處或高峰值函數(shù)附近小區(qū)域進(jìn)行粒子局部加密(見(jiàn)圖3(b)).

值得注意的是,在第4節(jié)和第5節(jié)的數(shù)值算例模擬中,均采用FPM離散(12)式結(jié)合3.2節(jié)邊界處理和時(shí)間項(xiàng)二階精度離散格式對(duì)C-H方程(1)—(4)式進(jìn)行離散求解.

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

本節(jié)數(shù)值求解有解析解但邊界條件不同的一/二維C-H方程,并結(jié)合數(shù)值結(jié)果分析LR-FPM方法的數(shù)值誤差和收斂速度.所模擬的問(wèn)題是一維帶周期邊值、二維帶Neumann邊值問(wèn)題,并與解析解做比較.為體現(xiàn)本文給出的方法相對(duì)于網(wǎng)格類(lèi)方法的優(yōu)點(diǎn),著重討論了在FPM方法的基礎(chǔ)上進(jìn)行局部加密分布和非均勻分布情況.局部加密采用了疏密點(diǎn)相結(jié)合的分布方式,這種處理方式的優(yōu)點(diǎn)在于: 1)相對(duì)于全局加密,局部加密增加數(shù)值精度的同時(shí)減少了計(jì)算量;2)相對(duì)于網(wǎng)格類(lèi)方法,局部加密簡(jiǎn)單易行.

為了分析LR-FPM方法的精度和收斂性,定義如下誤差(精確解與數(shù)值解L2-范數(shù)誤差)和收斂階:

其中 u 為數(shù)值解;U 為解析解;d1,d2為不同的粒子初始間距.

4.1 一維周期邊界C-H方程

考慮區(qū)域 Ω =[0,2π] 上帶周期邊值條件的C-H方程[7],其對(duì)應(yīng)的方程和初邊值條件分別為

該算例模擬中,均勻分布情況下粒子初始間距為 d1=π/32,時(shí)間步長(zhǎng)為 d t=10-4.局部加密情況下,在 [ 0 ,10π/32],[ 5 4π/32,2π] 處加密,加密區(qū)域的粒子間距為初始間距的1/2,即加密區(qū)域間距為d2=π/64,[ 1 0π/32,54π/32] 內(nèi)仍采用粒子初始間距d1,時(shí)間步長(zhǎng)為 d t=10-4,加密臨界點(diǎn)光滑長(zhǎng)度取h≈0.5×(1.1×d1+1.1×d2).圖 1展示了不同時(shí)刻下均勻分布和局部加密兩種情況下的FPM數(shù)值解,并與解析解作對(duì)比.可知兩種情況下的LRFPM結(jié)果均與解析解吻合.

圖1 幾個(gè)不同時(shí)刻下均勻分布、局部加密情況下的數(shù)值解和解析解Fig.1.Comparisons between the present numerical results and analytical solutions with different times under the uniform and local refinement particle distributions.

為體現(xiàn)提出的FPM方法模擬周期性問(wèn)題的數(shù)值精度和收斂速度,圖2給出了較短時(shí)間內(nèi)不同粒子數(shù)下數(shù)值誤差-時(shí)間曲線,L2-范數(shù)誤差 E2隨著粒子數(shù)的增加而減小,粒子數(shù)較多情況下具有更好的數(shù)值收斂性.值得注意的是在粒子數(shù)N=129時(shí)L2-范數(shù)誤差 E2先減小后增大,這是由空間步長(zhǎng)較小和隨時(shí)間延長(zhǎng)函數(shù)值變小情況下計(jì)算機(jī)的舍入誤差導(dǎo)致,但隨時(shí)間延長(zhǎng)中誤差 E2基本穩(wěn)定在同一量級(jí),仍符合誤差隨粒子增加而減小的數(shù)值計(jì)算理論.同時(shí)為了進(jìn)一步體現(xiàn)LR-FPM模擬周期性問(wèn)題的數(shù)值精度,表1列出了不同粒子間距下t=0.5s時(shí)的L2-范數(shù)誤差和收斂階.由表1可知,LR-FPM具有二階收斂速度,與文獻(xiàn)[7]中有限差分方法的收斂階基本一致,從而表明LR-FPM模擬C-H方程時(shí)具有較高的數(shù)值精度,較網(wǎng)格類(lèi)方法具有更好的靈活性和推廣應(yīng)用性.

圖2 不同粒子數(shù)下的收斂速度隨時(shí)間的變化Fig.2.The numerical convergence versus time under different particle numbers.

表1 t=0.5s 時(shí)不同粒子間距情況下的 L2-范數(shù)誤差 E2和收斂階Table 1.The L2-norm error E2 and convergence rate at t=0.5s .

表2列出了不同時(shí)刻均勻分布與局部加密情況的L2-范數(shù)誤差 E2,由表2可知局部加密情況下的誤差遠(yuǎn)遠(yuǎn)小于均勻分布情況下的誤差.這表明LR-FPM不僅具有更好的靈活推廣應(yīng)用性,而且能保持較高的精度.

表2 不同時(shí)刻下粒子均勻分布與局部加密情況下的L2-范數(shù)誤差 E2 對(duì)比Table 2.The L2-norm error E2 at different times under the uniform and local refinement particle distributions.

4.2 二維Neumann邊界C-H方程

為體現(xiàn)提出的LR-FPM方法求解帶第二類(lèi)邊界非線性C-H方程的準(zhǔn)確性,本節(jié)考慮正方形區(qū)域 Ω =[0,1]×[0,1] 的第二類(lèi)邊界條件的非線性C-H方程,其對(duì)應(yīng)的方程和初邊值條件[6]為

圖3展示了本算例模擬中采用的三種粒子分布方式,其中均勻分布方式是分別沿 x,y 方向每隔0.04的距離分布一個(gè)粒子,局部加密分布采用四角和中間加密,分別在 [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)為圓心,向外圓形分布粒子,相鄰的兩個(gè)圓距離相等,其余區(qū)域上仍采用均勻布點(diǎn)的方式.

圖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,時(shí)間步長(zhǎng)為 d t=10-6,局部加密采用加密區(qū)域粒子間距為 d2=0.025,其余區(qū)域粒子間距仍為 d1=0.04,加密臨界點(diǎn)光滑長(zhǎng)度取h≈0.5×(1.1×d1+1.1×d2).

圖4展示了 t=0.01s 時(shí)刻均勻分布和局部加密兩種粒子分布情況下 x=y 處的數(shù)值解,并與解析解進(jìn)行對(duì)比.可以看出,隨時(shí)間演化,兩種情況下的數(shù)值結(jié)果都與解析解吻合,但局部加密情況下的結(jié)果更接近解析解.

圖4 均勻分布與局部加密情況下的數(shù)值結(jié)果Fig.4.The present numerical results under the uniform and local refinement particle distributions.

表3列出了粒子間距為 d0=0.04 在FPM方法下采用五次樣條核函數(shù)與高斯核函數(shù)情況下的L2-范數(shù)誤差 E2對(duì)比.由表 3 可知,采用五次樣條核函數(shù)的誤差比采用高斯核函數(shù)的誤差小.因此,為提高數(shù)值精度,模擬中均采用五次樣條核函數(shù).

表3 初始間距 d0=0.04 情況下五次樣條核函數(shù)與高斯核函數(shù)的L2-范數(shù)誤差 E2 對(duì)比Table 3.The L2-norm error with quintic spline kernel and gaussian kernel functions at d0=0.04 .

為進(jìn)一步體現(xiàn)提出的LR-FPM方法模擬Neumann邊界問(wèn)題的數(shù)值精度和收斂速度,表4列出了模擬較短時(shí)間內(nèi)數(shù)值結(jié)果的誤差和收斂階,通過(guò)表4可得: 數(shù)值誤差隨著粒子數(shù)的增加而減小,且 d0=1/20 到 d0=1/40的誤差變化比d0=1/40 到 d0=1/60 的誤差變化大.給出的LRFPM方法模擬Neumann邊界條件下C-H方程具有二階收斂速度.

表5列出了粒子均勻分布、局部加密分布與非均勻分布情況下的L2-范數(shù)誤差 E2.通過(guò)觀察表5,本文給出的粒子方法在粒子分布均勻和非均勻情況下得到的數(shù)值誤差接近,局部加密情況下的誤差遠(yuǎn)遠(yuǎn)小于均勻分布情況下的誤差,表明局部加密能有效地提高數(shù)值精度.

表4 t=0.01s 時(shí)刻下不同粒子間距的 L2-范數(shù)誤差 E2和收斂階Table 4.The L2-norm error E2 and convergence rate at t=0.01s .

表5 粒子均勻分布、局部加密分布與非均勻分布情況下的L2-范數(shù)誤差 E2 對(duì)比Table 5.The L2-norm error E2 at different times under the uniform,local refinement,and non-uniform particle distributions.

為進(jìn)一步體現(xiàn)粒子方法在粒子非均勻分布情況下數(shù)值模擬精度,表6列出了 t=0.01s 時(shí)不同粒子間距非均勻分布情況下的誤差和收斂階.由表4、表5和表6可知,粒子方法在均勻分布和非均勻分布情況下得到的數(shù)值結(jié)果誤差都比較接近,表明了該方法易推廣應(yīng)用到非均勻區(qū)域問(wèn)題的模擬,且保持較好的精度.本文的粒子方法不受網(wǎng)格限制,可以在模擬區(qū)域任意粒子布置情況下對(duì)問(wèn)題進(jìn)行模擬實(shí)現(xiàn),較網(wǎng)格類(lèi)方法更容易被推廣應(yīng)用于復(fù)雜非規(guī)則區(qū)域上C-H問(wèn)題的模擬.

表6 t=0.01 s 時(shí)不同粒子間距非均勻分布情況下的L2-范數(shù)誤差 E2和收斂階Table 6.The L2-norm error E2 and convergence rate at t=0.01 s under non-uniform particle distribution.

5 C-H 方程數(shù)值模擬

為了更好地體現(xiàn)提出的LR-FPM粒子方法數(shù)值模擬C-H問(wèn)題的準(zhǔn)確性,本節(jié)分別選取一維/二維Neumann邊界無(wú)解析解C-H問(wèn)題,并與文獻(xiàn)[8]中基于離散變分導(dǎo)數(shù)的有限差分(finite difference method,FDM)法模擬結(jié)果做對(duì)比.

5.1 一維Neumann邊界C-H方程

考慮區(qū)域 Ω=[-1,1] 上帶有Neumann邊值條件的C-H方程,其表達(dá)式為

對(duì)應(yīng)的初值條件為u0(x)=0.53x+0.47sin(-1.5πx)和邊值條件為

該算例為一維無(wú)解析解的帶Neumann邊界C-H方程,它將描述復(fù)雜的相位分離現(xiàn)象,常被用來(lái)證明其保持方程的質(zhì)量守恒性質(zhì)和能量耗散性質(zhì).本節(jié)運(yùn)用LR-FPM方法對(duì) ε2=0.3 情況下算例進(jìn)行了模擬,并與文獻(xiàn)[8]中FDM方法結(jié)果進(jìn)行對(duì)比.圖 5 給出了 ε2=0.3,時(shí)間步長(zhǎng)為 d t=10-6,初始時(shí)刻與t=0.2s,t=2s 三個(gè)不同時(shí)刻兩種數(shù)值方法的模擬結(jié)果.由圖5可知,LR-FPM方法得到的帶Neumann邊界C-H方程的數(shù)值解與FDM 結(jié)果吻合.圖6給出了ε2=0.03,t=0.2s 時(shí)刻下粒子局部加密情況,在 [–1.0,–0.8],[0.8,1.0]處加密,每隔 0.02 的距離布置粒子,在(–0.8,0.8)內(nèi)每隔0.04的距離布置粒子.由圖6可知,粒子均勻分布與局部加密情況下的數(shù)值結(jié)果均與FDM結(jié)果吻合,但局部加密情況下的模擬結(jié)果更接近于FDM的結(jié)果.這進(jìn)一步表明提出的LR-FPM粒子方法模擬C-H方程是可靠的.

圖5 ε2=0.3 時(shí)不同時(shí)刻 FDM 結(jié)果與 LR-FPM 結(jié)果Fig.5.The numerical results obtained using FDM and LRFPM at different times with ε2=0.3 .

圖6 ε2=0.03,t=0.2s 時(shí)刻下均勻分布與局部加密情況下數(shù)值結(jié)果對(duì)比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方程

為體現(xiàn)提出的FPM方法求解帶Neumann邊界非線性C-H方程的準(zhǔn)確性,本節(jié)考慮正方形區(qū)域 Ω =[0,0.8]×[0,0.8] 的第二類(lèi)邊界條件的非線性C-H方程,其對(duì)應(yīng)的方程和初邊值條件[11]為

初值條件為

和Neumann邊界條件

圖7 t=0.1s 時(shí)刻 FPM 方法模擬結(jié)果Fig.7.The FPM result at t=0.1s .

圖7給出了 t=0.1s 時(shí)刻LR-FPM方法的模擬結(jié)果,與文獻(xiàn)[11]中數(shù)值結(jié)果相比較,可以看出兩種方法的結(jié)果接近.圖 8 給出了 t=0.08 s時(shí)刻GRBF方法[11]等值線圖(圖8(a))和本文方法在三種粒子分布情況下的數(shù)值等值線分布圖(圖8(b)—圖8(d)).通過(guò)觀察圖8可知,本文方法在粒子均勻分布、局部加密以及非均勻分布情況下得到的帶Neumann邊界C-H方程的數(shù)值等值線分布均與GRBF結(jié)果吻合,表明本文方法數(shù)值預(yù)測(cè)帶高階導(dǎo)數(shù)Neumann邊界C-H問(wèn)題是可靠的.

圖9給出了 t=0.08s 時(shí)刻下不同粒子間距均勻分布的情況以及局部加密分布下的數(shù)值收斂性.由圖9可知,LR-FPM數(shù)值預(yù)測(cè)二維C-H問(wèn)題時(shí)是收斂的.因此,LR-FPM能夠有效、可靠地模擬C-H方程,容易實(shí)施局部加密或非均勻分布情況下的數(shù)值模擬.

圖9 t=0.08s 時(shí)刻粒子局部加密分布情況下的數(shù)值收斂Fig.9.The numerical convergence obtained using the present method under different particle distributions at t=0.08s.

6 結(jié) 論

本文基于高階導(dǎo)數(shù)分解和區(qū)域離散局部加密思想,首次拓展應(yīng)用FPM方法對(duì)描述相位分離現(xiàn)象的含四階導(dǎo)數(shù)C-H方程進(jìn)行數(shù)值模擬研究,為復(fù)雜區(qū)域上含高階導(dǎo)數(shù)非線性問(wèn)題的數(shù)值預(yù)測(cè)提供了一種準(zhǔn)確、有效的LR-FPM方法.為了驗(yàn)證LR-FPM方法的優(yōu)勢(shì),首先采用LR-FPM求解了均勻、局部加密和非均勻離散情況下帶有解析解的一維/二維C-H方程,并分析了LR-FPM的數(shù)值誤差和收斂性;其次,運(yùn)用LR-FPM法對(duì)無(wú)解析解C-H方程進(jìn)行了數(shù)值預(yù)測(cè),并與有限差分法結(jié)果作對(duì)比.所有數(shù)值結(jié)果表明:

1)本文給出的帶有高階導(dǎo)數(shù)邊值條件的施加處理是準(zhǔn)確的,且給出的LR-FPM粒子方法模擬一維/二維C-H方程具有較好二階收斂速率;

2)基于局部加密下LR-FPM能有效地降低數(shù)值模擬誤差,且在粒子分布非均勻情況下比有限差分法有更好靈活應(yīng)用性;

3)給出的LR-FPM法能準(zhǔn)確預(yù)測(cè)不同物質(zhì)間相位分離非線性現(xiàn)象隨時(shí)間的演化過(guò)程.