羅 原 吳燕梅

(西南科技大學理學院 四川綿陽 621010)

晶體相場模型是由Elder等[1]提出的能夠在原子尺度和擴散時間上模擬材料微觀結(jié)構(gòu)的新方法, 是傳統(tǒng)的連續(xù)相場模型的延伸及推廣, 用以填補連續(xù)模型與分子動力學方法之間的空白。 該模型能夠耦合多重晶粒位向和彈塑性形變現(xiàn)象, 與傳統(tǒng)的密度泛函理論密切相關(guān)[2]，因此在研究晶體微觀結(jié)構(gòu)缺陷和大空間尺度下模擬材料微結(jié)構(gòu)變化更具優(yōu)勢。模型的能量泛函為：

(1)

其中，Ω=[0,L]×[0,L]∈R2,u:Ω→R表示密度, 一般情況下ε取小于1的常數(shù)。晶體相場方程的類型較多, 本文主要研究基于能量泛函式(1)在質(zhì)量守恒約束條件下導出的H-1梯度流：

ut=Δμ

(2)

(3)

假設u,Δu和μ在Ω上都滿足周期邊界條件，容易驗證能量泛函式(1)隨時間變化保持單調(diào)遞減。

由于晶體相場方程(2)是高階非線性偏微分方程, 如何突破對時間步長的苛刻限制是該方程數(shù)值模擬面臨的挑戰(zhàn)。Elder等[1]采用顯式歐拉格式進行模擬, 但顯式格式在穩(wěn)定性上要求時間步長Δt～O(h6)(h為空間步長)，導致此類格式對于長時間數(shù)值模擬基本是不可能的。全隱式格式雖然能增加時間步長, 但大部分算法在時間層上的精度不高。因此, 目前所研究的格式大部分都是半隱格式, 這樣既能夠在一定程度上放寬時間步長約束, 又能夠提高逼近精度, 同時又保證了數(shù)值格式的能量穩(wěn)定性。Backofen等[2]采用后向歐拉格式進行模擬,對非線性項采用線性化逼近，但沒有給出嚴格的穩(wěn)定性分析。Gomez等[3]利用修正的Crank-Nicolson方法，給出具有二階精度的非線性數(shù)值格式和相應的穩(wěn)定性分析。 Wise等[4]和Hu等[5]利用凸分解思想提出二階精度的有限差分格式，較完整證明了格式的質(zhì)量守恒性、數(shù)值解的存在唯一性以及格式的能量穩(wěn)定性和收斂性。Lee等[6]將晶體相場方程分解為線性和非線性兩部分,提出算子分裂格式,但在能量穩(wěn)定性上沒有給出具體分析。Vignal等[7]采用二階有限元方法, 并在理論上對格式的質(zhì)量守恒和能量穩(wěn)定性進行了較詳細分析。更多二階時間精度的數(shù)值算法的研究,可參看文獻[8-11]。為了較全面展現(xiàn)相變量的演化過程,晶體相場方程的數(shù)值模擬是相當耗費時間的,因此無論在時間上還是在空間上,高精度的數(shù)值格式對于方程的模擬是非常重要的。在時間高精度數(shù)值方法上, 相關(guān)的文獻并不多見。Shin等[12]利用凸分裂的思想給出了三階精度的凸分裂龍格-庫塔格式,同時分析了格式的能量穩(wěn)定性。

本文主要針對晶體相場方程(2)-方程(3)的離散，在時間上采用三階的后向差分公式，在空間上用Fourier 擬譜逼近，構(gòu)造時間和空間上都具有高階精度的數(shù)值格式。為保證計算格式的能量穩(wěn)定性,額外增加Douglas-Dupont正則項。在理論上對數(shù)值解的質(zhì)量守恒性、存在唯一性以及數(shù)值格式的能量穩(wěn)定性給出了詳細分析，數(shù)值算例驗證了格式的高精度以及晶體相場方程的長時間動力學模擬的有效性。

1 Fourier擬譜逼近

(4)

如果f∈GN, 它的有限離散Fourier展開式為：

(5)

(6)

(7)

顯然，▽N·▽Nf=ΔNf。

對網(wǎng)格函數(shù)f∈GN，引入離散算子(-ΔN)γ,γ>0，則:

(8)

(9)

利用離散的分部積分公式可以證明如下結(jié)論都是成立的：

(10)

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

(11)

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

(12)

基于上述定義, 可以得到晶體相場方程能量(1)的離散形式，即對于任何的周期網(wǎng)格函數(shù)u∈GN,離散能量為：

(13)

2 數(shù)值格式以及可解性和能量穩(wěn)定性

根據(jù)晶體相場模型的能量泛函結(jié)構(gòu), 基于凸分裂的思想, 把能量泛函分解為凹凸兩部分。晶體相場方程中對應的凸項采用隱式處理, 而凹項部分采用顯式處理。方程(2)-方程(3)中時間離散采用三階向后差分離散, 空間上用Fouier擬譜逼近。為保證數(shù)值格式的能量穩(wěn)定性和便于數(shù)值計算，對凹項顯式處理項進行外推, 增加Douglas-Dupont正則項, 得到如下具有精度O(Δt3+hm)的計算格式(m表示空間譜精度的階數(shù)):

(14)

(15)

(16)

式(14)-式(16)是一個三步格式, 而已知的初始層數(shù)據(jù)只有u0, 必須尋求其他的方法確定虛擬層u-1,u-2的值。如果設定u-1=u-2=u0, 能量穩(wěn)定性的分析比較容易，但初始層數(shù)據(jù)沒有誤差會導致計算結(jié)果在初始階段三階時間精度蛻化。當然, 也可以選擇其他高精度的格式如三階或者四階的龍格-庫塔格式計算u1,u2, 能夠保證精度不會損失。

證明：對于式(14)中的ΔNμn+1, 根據(jù)周期邊界條件, 利用離散分部積分公式可以得到：

(17)

定理2：如果數(shù)值格式式(14)- 式(16)中un,un-1,un-2∈GN, 則式(14)- 式(16)存在唯一解un+1。

證明：令

(18)

(19)

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

(20)

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

(21)

(22)

證明：將式(14)與(-ΔN)-1(un+1-un)作內(nèi)積, 得到：

(23)

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

(24)

采用離散的分部積分公式及等式(5), 對于非線性項, 得到估計：

女子處于土狼的包圍與攻擊中，跌跌撞撞，險象環(huán)生。體力的流逝，令出刀的速度和力量都大打折扣，甚至很難再對土狼造成致命的傷害。

(25)

類似式(16)的分析, 可以得到：

(27)

而且, 如下3個等式也是顯然成立的:

(28)

(29)

(30)

將式(25)- 式(30)代入式(23), 可得：

利用不等式：

(32)

3 數(shù)值算例

3.1 時間和空間精度測試

將空間二維計算區(qū)域設為Ω=[0,1]×[0,1], 相變量的精確解為:

(33)

此時方程(2)的右端不再是0，而是必須增加一個對應的外力項f(x,y,t), 即：

(34)

對時間層上的誤差估計, 需要把空間譜精度帶來的影響減少到最小。取空間方向網(wǎng)格劃分數(shù)N=256, 使得空間譜精度的誤差可以忽略。為比較不同參數(shù)組合下的誤差，選擇兩個參數(shù)A和ε的不同組合, 計算終止時間取t=1。為方便收斂階測試, 時間步長從0.1開始依次減半, 相變量u的數(shù)值解的‖·‖2誤差及收斂階數(shù)的計算結(jié)果見表1。

表1 不同參數(shù)下時間方向的‖·‖2誤差及階數(shù)Table 1 The‖·‖2 errors and orders in temporal direction with different parameters

從表1結(jié)果可以看出, 盡管理論上對系數(shù)A的要求極高, 在數(shù)值模擬時卻沒有這一限制條件。無論是A=1,A=10還是A=20,都沒有滿足定理3中的條件, 但是計算結(jié)果的精度非常高。另一方面, 隨著A值不斷增加, 絕對誤差越來越大, 這是因為A越大會產(chǎn)生更大的數(shù)值耗散。

對于空間方向的誤差估計, 選擇A=1,ε=0.025以及時間步長Δt=10-4。 時間步長取這么小的原因也是不希望時間層上的誤差給空間誤差的驗證帶來過多的干擾。 計算終止時間t=1，誤差結(jié)果見表2。顯然，F(xiàn)ourier譜精度隨著空間網(wǎng)格劃分數(shù)N的增加快速提升, 但是后期無論N值怎么增加, 精度都沒有大幅提高, 因為此時已經(jīng)不是空間離散誤差而是時間離散誤差占據(jù)主要作用。

表2 空間方向的‖·‖2誤差 Table 1 The ‖·‖2 errors for the spatial direction

3.2 晶體相場方程的動力學數(shù)值仿真

設二維空間區(qū)域為Ω=[0，128]×[0，128], 數(shù)值格式式(14)- 式(16)中的參數(shù)選取為A=2,ε=0.025, 空間步長固定為h=1, 晶體相的初始狀態(tài)值為：

ui,j=0.06+0.03ri,j

(35)

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

圖1 不同時刻晶體相場的狀態(tài)圖Fig.1 State diagram of the crystal phase field at different times

圖2給出不同時刻下的離散能量和離散質(zhì)量的計算結(jié)果。 可以看出離散能量的變化和理論一致, 即保持單調(diào)不增。前一段由于時間步長較小, 使得圖2中能量變化和時間步長較大的后一段相比較為平緩；中間段能量的變化較劇烈, 但隨著能量衰減到一定程度, 能量變化呈現(xiàn)出平緩狀, 此時晶體的相或者密度分布基本處于平穩(wěn)狀態(tài)。質(zhì)量隨著時間的變化，始終保持恒定，也驗證了定理1的正確性。

圖2 晶體離散能量和離散質(zhì)量的動態(tài)演化圖Fig.2 Dynamic evolution of the crystal discrete energy and discrete mass

4 結(jié)論

本文結(jié)合晶體相場模型的能量泛函結(jié)構(gòu), 利用凸分裂思想, 對周期邊界條件的晶體相場方程提出了一種具有能量穩(wěn)定性的數(shù)值算法。算法在空間上采用Fourier擬譜逼近，在時間采用上后向差分的三階精度逼近, 并增加Duglas-Dupont正則項以保證數(shù)值格式的能量穩(wěn)定性。對數(shù)值格式解的存在唯一性以及格式的能量穩(wěn)定性進行了證明，采用數(shù)值仿真驗證了格式在時間上的三階精度和空間上譜精度, 并給出了方程的長時間動力學演化模擬。