相場晶體方程的一個高精度能量穩(wěn)定數(shù)值格式

李貴川 , 賴 倩 , 胡勁松 , 馮端宇

(1.西華大學土木建筑與環(huán)境學院, 成都 610039；2. 西華大學理學院 成都, 610039； 3. 四川大學數(shù)學學院, 成都 610064)

1 引 言

作為一種在空間原子尺度上模擬晶體的新方法, 相場晶體模型[1]能夠通過極小化Swift-Hohenberg (SH) 型自由能量泛函來解釋晶格的周期性結(jié)構(gòu)及彈性和晶格的塑性變形、位錯、晶界等各種微結(jié)構(gòu)現(xiàn)象.此類模型主要采用相變量來描述原子密度的時間平均值,并與動態(tài)密度泛函理論密切相關(guān)[2]. 相場晶體模型具有能量泛函

(1)

其中Ω=[0,L]×[0,L]∈R2,φ:Ω→R表示密度場,a=1-ε為常數(shù),ε是與相場溫度有關(guān)的常數(shù). 相場晶體方程是由能量泛函(1)導(dǎo)出的H-1梯度流,即

φt=Δ(φ3+aφ+2Δφ+Δ2φ)

(2)

本文假設(shè)φ和Δφ在Ω上都滿足周期邊界條件.因此，方程(2)是質(zhì)量守恒的，同時其能量泛函(1)隨時間變化保持單調(diào)不增.

相場晶體方程(2)是一個六階非線性偏微分方程，難以獲得解析解，因而尋找高效數(shù)值算法尤為重要.構(gòu)造相場晶體方程數(shù)值算法的主要困難在于非線性項和高階導(dǎo)數(shù)項對時間步長存在苛刻的穩(wěn)定性限制.此外，由于需要在較長時間模擬過程中盡可能保證結(jié)果的準確性，數(shù)值格式的能量穩(wěn)定性也是一個關(guān)鍵因素.Backofen等[2]提出的有限元方法本質(zhì)上是標準后向歐拉格式，其中的非線性項φ3通過 (φk+1)3≈3(φk)2φk+1-2(φk)3逼近，因而無法保證能量穩(wěn)定性和可解性；Cheng和Warren[3]提出的線性化譜格式則不能在理論上證明該格式無條件能量穩(wěn)定；Yang等[4]給出的線性化一階及二階數(shù)值格式具有能量穩(wěn)定性，但φ的H2-穩(wěn)定性無法得到證明.此外，Hu及Wise等[5,6]利用凸分解方法討論了相場晶體模型的無條件穩(wěn)定的差分格式, 給出了時間方向一階及二階、空間方向二階精度的無條件收斂證明；李娟[7]提出了具有時間二階精度和空間四階精度的線性化三層差分格式，研究了格式的收斂性但沒有給出能量穩(wěn)定性分析. 另一方面，在時間高精度數(shù)值方法上,目前僅有的工作是Shin等[8]利用凸分解的思想構(gòu)造的時間上精度能達到三階的凸分裂龍格-庫塔格式，格式具有能量穩(wěn)定性和唯一可解性.更多的數(shù)值算法可參看文獻[9-12].

本文針對相場晶體方程在時間上采用三階后向差分公式、空間上利用理論上具有譜精度的Fourier擬譜逼近構(gòu)造了一個高精度數(shù)值格式.為保證數(shù)值格式的能量穩(wěn)定性，我們在格式中增加與時間階精度匹配的Douglas-Dupont正則項，進而在理論上證明了格式對于修正離散能量泛函的穩(wěn)定性.最后本文給出了數(shù)值算例，以驗證格式精度以及長時間數(shù)值模擬的有效性.

2 Fourier擬譜逼近及數(shù)值格式

為表述方便，本文選擇單位正方形區(qū)域，即L=1. 假設(shè)hx=hy=h,Nx=Ny=N,N=2K+1始終是奇數(shù). 網(wǎng)格點記為 (xi,yj),xi=ih,yj=jh, 0≤i,j≤N. 對于二維數(shù)值網(wǎng)格上的周期函數(shù)f, 定義網(wǎng)格函數(shù)空間

如果f,g∈GN，相應(yīng)的有限項離散Fourier展開式為

(3)

并且有?N·?Nf=ΔNf. 對網(wǎng)格函數(shù)f，引入離散算子(-ΔN)γ,γ>0，

對于定義的這些算子，通過離散的分部積分公式可以確保如下列等式恒成立：

〈f,ΔNg〉=-〈?Nf,?Ng〉，

(4)

〈f,g〉-1,N=〈f,(-ΔN)-1g〉=

〈(-ΔN)-1/2f,(-ΔN)-1/2g〉.

所以‖·‖-1,N可以定義為

至此，我們得到相場晶體方程能量的離散形式：

對于任何的周期網(wǎng)格函數(shù)φ∈GN，離散能量為

(5)

對于相場晶體方程(2)，時間離散采用三階向后差分離散，空間上采用Fouier擬譜逼近，同時為保證格式的能量穩(wěn)定性，增加Douglas-Dupont正則項，從而得到如下具有精度O(Δt3+hm)的數(shù)值格式:

(6)

此處

μn+1=(φn+1)3+aφn+1+2ΔN(3φn-3φn-1+

(7)

m表示空間譜精度的階數(shù)，具體取值依賴于空間網(wǎng)格點的劃分.

注1格式(6)是一個三步數(shù)值格式，初始層數(shù)據(jù)φ-1，φ-2必須給定. 如果選擇φ-1=φ-2=φ0, 能量估計會變得更簡單, 但是在初始階段，時間層的三階精度會退化. 另一方面, 也可以選擇三階的龍格-庫塔格式計算φ1,φ2，這樣就能夠保證初始時間層上達到理論上的三階收斂精度.

3 格式的可解性及能量穩(wěn)定性

證明 對于(6)式中的ΔNμn+1，根據(jù)周期邊界條件，可以得到

(8)

定理3.2如果φn,φn-1,φn-2∈GN，則數(shù)值格式(6)(7)存在唯一解φn+1.

證明 令

(9)

此處

(10)

非線性方程(9)的解能夠看作下述離散能量泛函的極小值：

根據(jù)凸優(yōu)化理論，該泛函是嚴格凸的，因此該泛函的極小值存在且唯一，也就表明格式(6)(7)的數(shù)值解也是存在唯一的. 證畢.

(11)

此處

(12)

證明 將式(6)與(-ΔN)-1(φn+1-φn)作內(nèi)積，得

〈ΔN(φn+1)3,(-ΔN)-1(φn+1-φn)〉-

a〈ΔNφn+1,(-ΔN)-1(φn+1-φn)〉-

2〈Δ2N(3φn-3φn-1+φn-2),(-ΔN)-1(φn+1-φn)〉-

AΔt2〈Δ2N(φn+1-φn),

(-ΔN)-1(φn+1-φn)〉=0

(13)

對上式中的第一項，可以推出

(14)

由離散的分部積分公式及式(4), 對于非線性項有

(15)

及

?N(φn+1-φn)〉=2〈?N(-φn+1+φn+1-2φn+φn-1-(φn-2φn-1+φn-2)),?N(φn+1-φn)〉≥

(16)

而且，如下三個等式也是顯然成立的:

(17)

(18)

(19)

將式(14)～(19)代入式(13), 可得

(20)

注意到

(21)

注2形式上定理3.3的結(jié)果表明本文只能夠證明數(shù)值格式(6)(7)對于修正的離散能量(12)而不是原始離散能量(5)的穩(wěn)定性. 實際上，如果假設(shè)φ-1=φ-2=φ0，則離散能量(12)的穩(wěn)定性能夠保證離散能量(5)的穩(wěn)定性. 這是因為

即任何時刻的能量都不會超過初始時刻的能量值.

4 數(shù)值算例

4.1 精度驗證

本節(jié)首先驗證數(shù)值格式(6)-(7)在時間方向上的三階精度.計算區(qū)域Ω=[0,1]×[0,1]，并假設(shè)相變量的精確解為:

(22)

為了與方程(2)保持一致，在方程(2)的右端需要增加一個外力項f(x,y,t)，即

f(x,y,t)=(φe)t-Δ(φe3+aφe+2Δφe+Δ2φe)

(23)

對于時間方向的誤差驗證，為了避免空間譜精度帶來的影響，取空間方向網(wǎng)格劃分數(shù)N=256，使得空間方向的誤差可以忽略不計. 在計算過程中，選擇兩個參數(shù)A和a的不同組合，計算結(jié)束時間為T=1，時間步長從0.1開始依次減半計算，相變量φe相應(yīng)的‖·‖2誤差及收斂階數(shù)的計算結(jié)果見表1.

表1 不同參數(shù)組合下的時間方向的‖·‖2誤差及階數(shù)

從表1結(jié)果可以看出，正如在注3中所說，盡管理論上對系數(shù)A的要求極高，但在數(shù)值模擬時卻沒有這一限制條件， 無論是A=1還是A=10. 不過需要提及的是, 從表1中的誤差結(jié)果看出，雖然A的值選擇越大越能夠在理論上保證能量穩(wěn)定性，但越大的A也會導(dǎo)致更大的數(shù)值耗散, 從而導(dǎo)致絕對誤差變大.

對于空間方向的誤差估計，選擇A=1,a=0.975以及時間步長Δt=10-4. 隨著空間網(wǎng)格逐漸加密，計算相變量在T=1時的誤差，計算結(jié)果見表2. 從結(jié)果可以看出，譜精度隨著空間網(wǎng)格劃分數(shù)N的增加快速提升，但無法得到類似時間方向的誤差的具體階數(shù). 這主要由于譜精度的作用, 隨著N的增大, 數(shù)值解的絕對誤差已經(jīng)變化不大，因為這個時候的誤差已經(jīng)不是空間離散主要引起的，而是時間方向誤差的體現(xiàn).

表2 空間方向的‖·‖2誤差

4.2 相場晶體的動力學模擬

考慮二維空間區(qū)域Ω=[0,128]×[0,128], 數(shù)值格式(6)(7)中的參數(shù)選取為A=1,a=0.975，空間步長固定為h=1，晶體相的初始狀態(tài)值為

φi,j=0.06+0.03ri,j

(24)

其中ri,j為[0,1]服從均勻分布的隨機數(shù). 為提高計算效率，采取不同的時間步長進行數(shù)值模擬，即在時間區(qū)間[0,100]上取Δt=0.01，在[100,2000]上取Δt=0.1，相場晶體的動態(tài)演化在不同時間的狀態(tài)見圖1.

(a)

(b)

(c)

(d)

圖2也給出不同時刻下的離散能量泛函(5)的計算結(jié)果. 顯然，離散能量的演化和理論分析結(jié)果一致，即保持單調(diào)不增. 開始時由于時間步長較小，使得圖2中能量變化比后一段較大時間步長時能量變化平緩；中間段能量變化較劇烈，但隨著能量衰減到一定程度，能量變化呈現(xiàn)出平緩狀，此時晶體相或者密度分布基本處于平穩(wěn)狀態(tài).

接下來取晶體相的另外一種初始狀態(tài)下晶體微結(jié)構(gòu)的變化. 先取晶體相的初始狀態(tài)值為φi,j=0.02ri,j，然后再四個空間點處增加四個點狀相場值(見圖3a)，其值分別為10、5、8和10. 隨著時間增加，晶體相場的演化見圖3. 可以看出，在初期階段時四個點源項在演化過程中互不干擾，但較長時候后，各點的演化互相影響，同時由于相場值的差異，導(dǎo)致不對稱的微結(jié)構(gòu)出現(xiàn).

圖2 晶體能量的動態(tài)演化圖Fig.2 The evolution of the discrete energy

(a)

(b)

(c)

(d)

5 結(jié) 論

本文針對相場晶體方程的周期邊界問題，提出了一種具有能量穩(wěn)定性的Fourier擬譜數(shù)值格式. 由于能量穩(wěn)定性是在材料微觀結(jié)構(gòu)的長時間數(shù)值模擬時本質(zhì)物理屬性的要求, 因此我們在時間采用后向差分的三階高精度離散, 并增加Duglas-Dupont正則項以確保能量穩(wěn)定，并以理論上證明了數(shù)值格式的質(zhì)量守恒性、數(shù)值解的存在唯一性以及格式的能量穩(wěn)定性. 數(shù)值模擬結(jié)果既驗證了本格式在時間和空間上高精度離散逼近, 同時也展現(xiàn)了相場晶體方程所描述的材料結(jié)構(gòu)的微觀變化以及能量的演化趨勢.