羅 原 吳燕梅

(西南科技大學(xué)理學(xué)院 四川綿陽(yáng) 621010)

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

(1)

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

ut=Δμ

(2)

(3)

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

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

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

1 Fourier擬譜逼近

(4)

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

(5)

(6)

(7)

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

對(duì)網(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)

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

(13)

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

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

(14)

(15)

(16)

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

證明：對(duì)于式(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)化理論, 該泛函是嚴(yán)格凸的, 因此該泛函的極小值是存在且唯一的, 也就表明式(14)- 式(16)的數(shù)值解也是存在且唯一的。證畢。

(21)

(22)

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

(23)

對(duì)上式中的第一項(xiàng), 可以推出：

(24)

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

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

(25)

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

(27)

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

(28)

(29)

(30)

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

利用不等式：

(32)

3 數(shù)值算例

3.1 時(shí)間和空間精度測(cè)試

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

(33)

此時(shí)方程(2)的右端不再是0，而是必須增加一個(gè)對(duì)應(yīng)的外力項(xiàng)f(x,y,t), 即：

(34)

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

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

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

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

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

3.2 晶體相場(chǎng)方程的動(dòng)力學(xué)數(shù)值仿真

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

ui,j=0.06+0.03ri,j

(35)

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

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

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

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

4 結(jié)論

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