孫竟巍 張弘 錢旭 宋松和

(國防科技大學(xué)文理學(xué)院,長沙 410073)

1 引言

經(jīng)典Allen–Cahn 方程是相場模型中關(guān)于自由能泛函:

的一個典型的梯度流問題.這里,u為一個實值未知函數(shù),? 為具有Lipschitz 邊界?? 的有界區(qū)域.以(·,·)表示? 上的L2內(nèi)積,對應(yīng)的L2范數(shù)記為∥·∥2.

在兩相流問題中,u通常表示相變量,?>0 代表兩相間的界面寬度,F(u)則表示某種非線性勢函數(shù).經(jīng)典Allen–Cahn 方程可以視為在能量泛函(1.1)下的L2梯度流:

其中f(u)=?F′(u).方程(1.2)是由Samuel M.Allen 和John W.Cahn[1]于1979 年提出,用于描述結(jié)晶體中反相邊界的運動模型,如今已被廣泛用于模擬自然界中的界面運動和相變,如[2,3].Allen–Cahn 方程(1.2)具有能量耗散律:

此外,Allen–Cahn 方程還有一個重要的性質(zhì):極值原理(Maximum principle preserving,MPP).極值原理作為研究相場模型物理特性的一個重要的數(shù)學(xué)工具,在近些年受到越來越多學(xué)者的關(guān)注.極值原理是指如果方程的初值和邊值以某個大于0 的特定常數(shù)逐點為界,那么計算所得到的解也始終以該常數(shù)為界.作為相場模型的一個重要的物理特征,極值原理對于具有物理意義的數(shù)值解而言是必不可少的.不保持極值原理的數(shù)值算法可能使得數(shù)值解出現(xiàn)病態(tài)情況并且數(shù)值模擬會發(fā)生“爆破”現(xiàn)象.所以人們在對Allen–Cahn 方程所描述的動力學(xué)進行數(shù)值模擬時,希望保持這個特性.

經(jīng)典Allen–Cahn 方程雖然滿足能量耗散和極值原理,但由于其并不遵循質(zhì)量守恒定律,可研究的范圍具有一定的局限性.為進一步改善Allen–Cahn 模型來擴展相關(guān)動力學(xué)研究范圍,學(xué)者進行了廣泛且深入的研究.目前主要有兩類修正的守恒型Allen–Cahn 方程.一類是由Rubinstein 和Sternberg[4]將非局部具有空間依賴的拉格朗日算子引入到經(jīng)典Allen–Cahn 方程中,使得修正的方程解具有能量守恒性質(zhì)?該算子簡記為RSLM,形式為:

另一類是由Brassel 和Bretin[5]將局部-非局部且時空依賴的拉格朗日乘子引入到經(jīng)典Allen–Cahn方程使其解滿足質(zhì)量守恒?該算子簡記為BBLM,形式為:

近幾年內(nèi),學(xué)者對Allen–Cahn 型方程極值原理這個特性的關(guān)注密切,許多針對該方程構(gòu)造的保極值格式應(yīng)運而生.文獻[7,8]針對經(jīng)典Allen–Cahn 方程將有限差分法與穩(wěn)定化半隱方法相結(jié)合構(gòu)造了一階和二階的保極值格式.杜強等[9]針對非局部Allen–Cahn 方程構(gòu)造了一階和二階的指數(shù)時間差分格式.侯天亮等[10] 對空間分?jǐn)?shù)階Allen–Cahn 方程提出了保極值的二階Crank–Nicolson 有限差分方法.最近,廖洪林等[11]針對時間分?jǐn)?shù)階Allen–Cahn 方程提出了二階的時間自適應(yīng)保極值格式.然而很少有文獻研究關(guān)于守恒型Allen–Cahn 方程的保極值格式.文獻[12]針對守恒型Allen–Cahn 方程構(gòu)造了無條件保極值格式,但該格式精度最高只能達到二階.就我們所知,還沒有任何文獻針對守恒型Allen–Cahn 方程構(gòu)造的保極值格式精度超過二階.而高階精度對于高維和強復(fù)雜度問題的研究十分重要,因此我們迫切需要一個高精度的保極值格式來求解守恒型Allen–Cahn 方程.

本文借助時間方向上的積分因子兩步Runge–Kutta 方法來構(gòu)造高精度保極值格式.積分因子Runge–Kutta 方法近些年來已被廣泛應(yīng)用到剛性常微分方程的計算中[13,14],積分因子Runge–Kutta 方法可以看作是傳統(tǒng)Runge–Kutta 方法的一個延伸,尤其對于含有強剛性的線性項問題格外有效.其主要思想在于應(yīng)用指數(shù)積分因子來消除線性項的強剛性并且用Runge–Kutta 方法進行時間積分.文獻[15,16]使用積分因子Runge–Kutta 方法發(fā)展了經(jīng)典Allen–Cahn 方程的最高四階精度保極值格式.文獻[17]針對經(jīng)典Allen–Cahn 方程將積分因子Runge–Kutta 方法結(jié)合線性穩(wěn)定化技術(shù)構(gòu)造了最高三階精度的無條件保極值格式.本文使用的兩步Runge–Kutta 方法(Two-step Runge–Kutta method,TSRK)作為傳統(tǒng)Runge–Kutta 方法的擴展,有接近30 年的研究歷史.兩個方法的區(qū)別在于兩步Runge–Kutta 方法應(yīng)用前一時間層和當(dāng)前時間層的信息構(gòu)造積分器,繼而使得兩步Runge–Kutta 法可以用較少的級數(shù)得到與單步Runge–Kutta 法相同的精度.

本文的其余部分安排如下.第二節(jié)介紹關(guān)于RSLM 形式的守恒型Allen–Cahn 方程和半離散系統(tǒng)的一些預(yù)備知識.第三節(jié)構(gòu)造積分因子兩步Runge–Kutta 格式并證明該格式可以保持極值原理和保持質(zhì)量守恒,隨后給出該格式的誤差分析.在第四節(jié),通過一系列的數(shù)值實驗來驗證數(shù)值格式的有效性和優(yōu)勢所在.在最后給出結(jié)論和展望.

2 預(yù)備知識

2.1 具有非局部拉格朗日乘子的守恒型Allen–Cahn 方程

在周期邊界條件下,將RSLM 形式的拉格朗日乘子[4]引入經(jīng)典Allen–Cahn 方程(1.2)中,得到RSLM 形式的守恒型Allen–Cahn 方程:

這里,修正的非線性項定義為

其中λ(u,t)=(|?|是? 的Lebesgue 測度)是拉格朗日乘子.對于RSLM 形式的守恒型Allen–Cahn 方程(2.1)而言,文獻[4,12]推導(dǎo)出(2.1)具有極值原理時對非線性項的要求.

假設(shè)1([4,12]) 存在一個常數(shù)β >0 使得

下面給出守恒型Allen–Cahn 方程(2.1)的極值原理.

定理1([4,12]) 給定一個常數(shù)T >0,假設(shè)u(t)∈C1(0,T;C2(?))是帶有周期邊界條件的守恒型Allen–Cahn 方程(2.1) 的經(jīng)典解.如果假設(shè)1成立并且對于任意x ∈? ∪??,初值滿足≤β,那么對于任意t>0,≤β成立.

對于廣泛使用的多項式勢函數(shù)F(u)=,非線性項為

文獻[12]證明了對于任意β ≥,(2.4)中的f(u)都滿足假設(shè)1.而當(dāng)β=時,我們有|f′(ξ)| ≤3,?|ξ| ≤β.因此由于非局部RSLM 算子λ(u,t)的引入不能保證解的局部性,從而改變了經(jīng)典Allen–Cahn 方程原本的極值.然而,在周期邊界條件下,這個非局部項抵消了質(zhì)量的變化,即

并且沒有影響原始的能量耗散律:

這里,能量泛函是(1.1),g(u)=

2.2 有限差分半離散系統(tǒng)

作為一個常用的空間離散方法[15,16,17],二階有限差分法將拉普拉斯算子離散為M-矩陣,這樣可以繼承格式中一些諸如極值原理等的關(guān)鍵性質(zhì).下面給出一維空間中二階有限差分離散的一些重要結(jié)果.

在周期邊界條件下考慮計算區(qū)間[xL,xR],其中網(wǎng)格尺寸h=網(wǎng)格點為?h={xj|xj=xL+jh,j=0,1,···,N ?1}.記VN={v|v=(vj),vj ∈?h} ?RN為定義在?h上的網(wǎng)格函數(shù)空間,其離散內(nèi)積和范數(shù)定義如下:

二階拉普拉斯算子?xx的中心差分矩陣記為

下面給出中心差分矩陣D2的一些基本結(jié)果.

引理1([18]) 關(guān)于矩陣D2,存在這樣一種關(guān)系

令L=?2D2.根據(jù)[19]的假設(shè)2,由離散矩陣L可以得到RN上的一個收縮半群,滿足

其中∥·∥∞為矩陣的無窮范數(shù).

注意到D2矩陣的所有行和都為0,我們有下面這條引理.

引理2([20]) 對于任意τ >0,總有=1.

引理3([20]) 對于矩陣L,有:

這里1=[1,1,···,1]T ∈RN.

一維空間下守恒型的Allen–Cahn 方程(2.1)可以寫為如下形式:

運用文獻[19]中的抽象框架,半離散系統(tǒng)(2.10)滿足下面的極值原理:

定理2([19]) 考慮勢函數(shù)為(2.4)的半離散系統(tǒng)(2.10).如果初值u0=u0(x)滿足≤β,那么(2.10)的數(shù)值解u(t)滿足離散極值原理

將離散質(zhì)量和能量分別記為Mh(u)和Eh(u).分別對(2.10)的兩邊與1 和ut作l2內(nèi)積并使用周期邊界條件,可得到半離散系統(tǒng)(2.10)的質(zhì)量守恒律:

和能量耗散律:

3 時間離散方法

本節(jié)結(jié)合時間方向上的積分因子兩步Runge–Kutta 方法對半離散系統(tǒng)(2.10)構(gòu)造全離散格式,并證明該格式可以保持質(zhì)量守恒和極值原理.最后給出該格式的誤差分析.

3.1 保極值的積分因子兩步Runge–Kutta 格式

保持強穩(wěn)定性(Strong stability-preserving,SSP) 的格式[21,22] 已經(jīng)被廣泛用來構(gòu)造保極值的數(shù)值格式.本文將這種格式與積分因子兩步Runge–Kutta 方法相結(jié)合,所得到的保持強穩(wěn)定性積分因子兩步Runge–Kutta(Strong stability-preserving integrating factor two-step Runge–Kutta,SSP-IFTSRK)格式可以在某種條件下保持(2.10)的極值原理和質(zhì)量守恒.

Isherwood 等在[23]中開發(fā)并分析了一系列最高到八階的帶有優(yōu)化SSP 系數(shù)的IFTSRK 格式.下面考慮一個帶有SSP 系數(shù)C的s級p階顯式TSRK 格式,定義其Butcher 表如下:

其中ai,j=0(i ≤j),ci=,i=?1,0,···,s并且所有的系數(shù)都受精度和穩(wěn)定性的限制.

在方程(2.10)的左右兩邊乘以e?Lt,可得到

因此,

應(yīng)用變換w=e?Ltu到上面的方程,(2.10)化為:

將(3.4)與TSRK 結(jié)合,可得到如下的IFTSRK 格式:

再將(3.5)重寫為Shu–Osher 的形式[23]:

運用定理2,可以推導(dǎo)出下面的定理.

定理3假設(shè)使用到的TSRK 格式有非負的系數(shù)非減的節(jié)點ci和SSP 系數(shù)C.在條件τ ≤下,其中τF E表示滿足

證明用數(shù)學(xué)歸納法來證明該定理.當(dāng)n=0 時,結(jié)論是顯然的.假設(shè)≤β,我們將驗證對于un+1結(jié)果依然成立.

對于(3.6)的每一級,有

對于任何τ ≤都成立.定理得證.

下面給出SSP-IFTSRK 格式(3.5)保持半離散系統(tǒng)(2.10)質(zhì)量守恒的定理及證明.

定理4假設(shè)TSRK 格式有非負的系數(shù)di,j和ai,j,非減的節(jié)點ci和SSP 系數(shù)C,那么在條件τ ≤CτF E下,SSP-IFTSRK 格式(3.5)可以保持守恒型Allen–Cahn 方程(2.1)的數(shù)值解的質(zhì)量守恒律,即如果那么?n>0.

證明同樣用數(shù)學(xué)歸納法來證明定理的結(jié)論.對于守恒型Allen–Cahn 方程(2.1),有下式成立:

由〈un+1,1〉=〈un,s,1〉可知定理結(jié)論成立.

3.2 誤差分析

本小節(jié)將建立關(guān)于守恒型Allen–Cahn 方程(2.1)的SSP-IFTSRK 格式在無窮范數(shù)意義下的數(shù)值解誤差分析.為了簡便起見,省略TSRK 方法的相關(guān)階條件.通過使用定理3和s級p階TSRK格式的定義,可以推導(dǎo)出如下的收斂結(jié)果.

定理5給定T >0,假設(shè)u(t)∈Cp[0,T]是半離散系統(tǒng)(2.10)的精確解,un是由SSP-IFTSRK方法(3.5)計算所得的數(shù)值解并且初值在[0,T]上充分光滑.在定理3的條件下,如果那么有誤差估計

證明參考[15,24] 的證明方法,假設(shè)參考解為Un,i,?1≤i ≤s,其中,Un,?1=u(tn?1),Un,0=u(tn),Un,s=u(tn+1),且

截斷誤差Rs滿足

其中Cs與u的Cp+1[0,T]范數(shù)及f,s,p,T的Cp[?β,β]范數(shù)有關(guān),但與τ無關(guān).

令en+j=u(tn+j)?un,j,j=?1,0,en,i=Un,i ?un,i,?1≤i ≤s,則en,j=en+j,j=?1,0和en,s=en+1.對于每個i=1,2,···,s ?1,由(3.9)和(3.5)作差可得

類似地,有

遞歸可得

4 數(shù)值實驗

本節(jié)首先介紹一些必要的準(zhǔn)備工作,之后測試格式(3.5)的時間收斂階,并且比較同一階數(shù)下不同級數(shù)的情況.最后分別在二維和三維情形下模擬守恒型Allen–Cahn 方程的相變現(xiàn)象并且測試能量耗散性,保極值特性和質(zhì)量守恒律.

4.1 準(zhǔn)備工作

為便于之后的數(shù)值測試,將幾種IFTSRK 方法和對應(yīng)的SSP 系數(shù)列在表1 中.

表1 幾個帶有非負系數(shù)的優(yōu)化SSP-IFTSRK(s,p)方法及其SSP 系數(shù)[23]

考慮帶有周期邊界條件的方程(2.1)時,d 維的DFT 技術(shù)可以應(yīng)用到指數(shù)算子上.對任意的u ∈RN,算子可以重寫為

其中F和F?1分別表示快速傅里葉變換(Fast Fourier transform,FFT)和它的逆變換,Λ由式(2.8)給出.

因為TSRK 方法不能自啟動而且涉及到兩個時間層上的變量,因此在用格式(3.5)進行數(shù)值測試前需要使用一個啟動格式.該啟動格式應(yīng)滿足兩個條件:具有足夠的精度和保極值特性.將第一個時間步一分為二后,使用SSPRK(5,4)格式或者SSPRK(10,4)格式作為第一個子時間步的計算,而第二個子時間步則用TSRK 格式進行計算.在之后的時間步里,只需用到時間層tn和tn+1對應(yīng)的值un?1和un即可.

4.2 收斂階測試

例1考慮以下一維初值條件[16]:

設(shè)參數(shù)?=0.01,h=下面,計算時間步長為τ=2?k,k=0,1,···,8 的數(shù)值解.參考解則由時間步長取τref=2?12的SSP-IFTSRK(11,8)格式計算得出.圖1的黑色虛線作為基準(zhǔn)線用來表示理論階的參考.

圖1 T=2 時不同數(shù)值格式間的精度比較

圖1 的左邊展示出SSP-IFTSRK 格式從五階到八階的數(shù)值誤差.四條線大致與基準(zhǔn)線保持平行.圖1 的右邊則展示出同一階數(shù)對應(yīng)不同級數(shù)的收斂情況.實驗與理論結(jié)果保持一致,即隨著階數(shù)增加,極值誤差逐漸減少.

4.3 保結(jié)構(gòu)測試

例2考慮常見的粘連氣泡標(biāo)準(zhǔn)問題[25],其初值條件設(shè)置為

另外,半徑R設(shè)置為0.19,(x1,y1),(x2,y2)分別設(shè)置為(0,0.2)和(0,?0.2),界面參數(shù)?=0.01,時間步長τ=0.5,網(wǎng)格點數(shù)NX=NY=128.

使用經(jīng)典AC 方程(1.2)和帶有RSLM 形式的守恒型方程(2.1)分別計算該問題.選取級數(shù)較少的SSP-IFTSRK(4,5)格式進行測試.另外,選取6 個時間節(jié)點并在圖2 展示出數(shù)值解的演化過程.圖3 則展示對兩個方程進行數(shù)值計算后得到的質(zhì)量、極值和能量變化情況.在守恒型Allen–Cahn方程的極值圖中我們添加一條黃色實線作為理論極值基準(zhǔn)線,值為

從圖2 可以看出,隨著時間的推移,經(jīng)典Allen–Cahn 方程對應(yīng)的兩個氣泡從粘連到結(jié)合,變成橢圓后形狀逐漸變小,最后在T=400 時徹底消失.而守恒型Allen–Cahn 方程對應(yīng)的兩個氣泡在結(jié)合后并沒有逐漸變小,形狀卻變得較為飽滿,最后處于一個穩(wěn)定的形態(tài).這說明用SSP-IFTSRK格式可以很好地模擬出傳質(zhì)現(xiàn)象并且能直觀地展示出兩個不同方程間的差異.

從圖3 可以看出,對于守恒型Allen–Cahn 方程,SSP-IFTSRK(4,5)格式能很好地保持質(zhì)量守恒,同時極值也始終在理論值內(nèi),能量曲線單調(diào)下降.而對于經(jīng)典Allen–Cahn 方程,可以看到質(zhì)量曲線單調(diào)下降,極值與理論值1 保持一致,能量曲線單調(diào)下降且在最后時刻降為0.三個量的變化情況與圖2 的演化過程是完全對應(yīng)的,這也體現(xiàn)出所提出的數(shù)值格式的有效性和魯棒性.

圖2 使用SSP-IFTSRK(4,5)計算的粘連氣泡數(shù)值解的演化過程.上半部分:經(jīng)典Allen–Cahn 方程?下半部分:守恒型Allen–Cahn 方程

圖3 使用SSP-IFTSRK(4,5)計算的粘連氣泡的質(zhì)量、極值和能量變化情況.上半部分:經(jīng)典Allen–Cahn 方程?下半部分:守恒型Allen–Cahn 方程

例3考慮隨機型三維初值問題:

設(shè)置界面參數(shù)?=0.01,網(wǎng)格點數(shù)NX=NY=NZ=128,時間步長τ=0.5.依然采用SSP-IFTSRK(4,5)格式對守恒型Allen-Cahn 方程(2.1)進行數(shù)值測試.選取6 個時刻來刻畫兩相分離的過程,分別為T=0,5,10,50,100,200.

數(shù)值解的演化過程在圖4 中展示.可以看到,隨著時間的推移,細小雜亂無章的相逐漸聚集并粗化,在T=200 時相的演化與初始時刻相比已經(jīng)截然不同.

圖4 使用SSP-IFTSRK(4,5)格式計算的三維隨機型數(shù)值解的演化過程

整個過程的質(zhì)量,極值和能量變化情況則在圖5 中展示,其中黃色實線表示理論值可以看到質(zhì)量是守恒的,能量曲線隨著時間單調(diào)下降,而且極值原理也保持得非常理想.

圖5 使用SSP-IFTSRK(4,5)格式計算的三維隨機型數(shù)值解的質(zhì)量,極值和能量變化過程

5 結(jié)論

本文構(gòu)造并分析了針對RSLM 形式的守恒型Allen–Cahn 方程的高階保極值數(shù)值格式,并且分別從理論和實驗兩個方面證明了所提出的數(shù)值格式可以保持該方程的質(zhì)量守恒,極值原理和能量耗散.數(shù)值結(jié)果也展現(xiàn)出所提出格式的一些優(yōu)勢,例如長時間穩(wěn)定性、高效性和魯棒性等.

注意到在使用SSP-IFTSRK 格式進行計算時,時間步長需滿足不大于CτF E這一條件.因此如何消除這個條件使其能夠無條件保極值原理,并且保持守恒型Allen–Cahn 方程的其他物理特性則是接下來要研究的目標(biāo).