晶体相场方程的高精度 Fourier 拟谱逼近

罗 原 吴燕梅

(西南科技大学理学院 四川绵阳 621010)

晶体相场模型是由Elder等[1]提出的能够在原子尺度和扩散时间上模拟材料微观结构的新方法, 是传统的连续相场模型的延伸及推广, 用以填补连续模型与分子动力学方法之间的空白。 该模型能够耦合多重晶粒位向和弹塑性形变现象, 与传统的密度泛函理论密切相关[2]，因此在研究晶体微观结构缺陷和大空间尺度下模拟材料微结构变化更具优势。模型的能量泛函为：

(1)

其中，Ω=[0,L]×[0,L]∈R2,u:Ω→R表示密度, 一般情况下ε取小于1的常数。晶体相场方程的类型较多, 本文主要研究基于能量泛函式(1)在质量守恒约束条件下导出的H-1梯度流：

ut=Δμ

(2)

(3)

假设u,Δu和μ在Ω上都满足周期边界条件，容易验证能量泛函式(1)随时间变化保持单调递减。

由于晶体相场方程(2)是高阶非线性偏微分方程, 如何突破对时间步长的苛刻限制是该方程数值模拟面临的挑战。Elder等[1]采用显式欧拉格式进行模拟, 但显式格式在稳定性上要求时间步长Δt～O(h6)(h为空间步长)，导致此类格式对于长时间数值模拟基本是不可能的。全隐式格式虽然能增加时间步长, 但大部分算法在时间层上的精度不高。因此, 目前所研究的格式大部分都是半隐格式, 这样既能够在一定程度上放宽时间步长约束, 又能够提高逼近精度, 同时又保证了数值格式的能量稳定性。Backofen等[2]采用后向欧拉格式进行模拟,对非线性项采用线性化逼近，但没有给出严格的稳定性分析。Gomez等[3]利用修正的Crank-Nicolson方法，给出具有二阶精度的非线性数值格式和相应的稳定性分析。 Wise等[4]和Hu等[5]利用凸分解思想提出二阶精度的有限差分格式，较完整证明了格式的质量守恒性、数值解的存在唯一性以及格式的能量稳定性和收敛性。Lee等[6]将晶体相场方程分解为线性和非线性两部分,提出算子分裂格式,但在能量稳定性上没有给出具体分析。Vignal等[7]采用二阶有限元方法, 并在理论上对格式的质量守恒和能量稳定性进行了较详细分析。更多二阶时间精度的数值算法的研究,可参看文献[8-11]。为了较全面展现相变量的演化过程,晶体相场方程的数值模拟是相当耗费时间的,因此无论在时间上还是在空间上,高精度的数值格式对于方程的模拟是非常重要的。在时间高精度数值方法上, 相关的文献并不多见。Shin等[12]利用凸分裂的思想给出了三阶精度的凸分裂龙格-库塔格式,同时分析了格式的能量稳定性。

本文主要针对晶体相场方程(2)-方程(3)的离散，在时间上采用三阶的后向差分公式，在空间上用Fourier 拟谱逼近，构造时间和空间上都具有高阶精度的数值格式。为保证计算格式的能量稳定性,额外增加Douglas-Dupont正则项。在理论上对数值解的质量守恒性、存在唯一性以及数值格式的能量稳定性给出了详细分析，数值算例验证了格式的高精度以及晶体相场方程的长时间动力学模拟的有效性。

1 Fourier拟谱逼近

(4)

如果f∈GN, 它的有限离散Fourier展开式为：

(5)

(6)

(7)

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

对网格函数f∈GN，引入离散算子(-ΔN)γ,γ>0，则:

(8)

(9)

利用离散的分部积分公式可以证明如下结论都是成立的：

(10)

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

(11)

所以‖·‖-1,N可以定义为：

(12)

基于上述定义, 可以得到晶体相场方程能量(1)的离散形式，即对于任何的周期网格函数u∈GN,离散能量为：

(13)

2 数值格式以及可解性和能量稳定性

根据晶体相场模型的能量泛函结构, 基于凸分裂的思想, 把能量泛函分解为凹凸两部分。晶体相场方程中对应的凸项采用隐式处理, 而凹项部分采用显式处理。方程(2)-方程(3)中时间离散采用三阶向后差分离散, 空间上用Fouier拟谱逼近。为保证数值格式的能量稳定性和便于数值计算，对凹项显式处理项进行外推, 增加Douglas-Dupont正则项, 得到如下具有精度O(Δt3+hm)的计算格式(m表示空间谱精度的阶数):

(14)

(15)

(16)

式(14)-式(16)是一个三步格式, 而已知的初始层数据只有u0, 必须寻求其他的方法确定虚拟层u-1,u-2的值。如果设定u-1=u-2=u0, 能量稳定性的分析比较容易，但初始层数据没有误差会导致计算结果在初始阶段三阶时间精度蜕化。当然, 也可以选择其他高精度的格式如三阶或者四阶的龙格-库塔格式计算u1,u2, 能够保证精度不会损失。

证明：对于式(14)中的ΔNμn+1, 根据周期边界条件, 利用离散分部积分公式可以得到：

(17)

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

证明：令

(18)

(19)

非线性方程(18)能够看作下述离散能量泛函的极小值：

(20)

根据凸优化理论, 该泛函是严格凸的, 因此该泛函的极小值是存在且唯一的, 也就表明式(14)- 式(16)的数值解也是存在且唯一的。证毕。

(21)

(22)

证明：将式(14)与(-ΔN)-1(un+1-un)作内积, 得到：

(23)

对上式中的第一项, 可以推出：

(24)

采用离散的分部积分公式及等式(5), 对于非线性项, 得到估计：

女子处于土狼的包围与攻击中，跌跌撞撞，险象环生。体力的流逝，令出刀的速度和力量都大打折扣，甚至很难再对土狼造成致命的伤害。

(25)

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

(27)

而且, 如下3个等式也是显然成立的:

(28)

(29)

(30)

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

利用不等式：

(32)

3 数值算例

3.1 时间和空间精度测试

将空间二维计算区域设为Ω=[0,1]×[0,1], 相变量的精确解为:

(33)

此时方程(2)的右端不再是0，而是必须增加一个对应的外力项f(x,y,t), 即：

(34)

对时间层上的误差估计, 需要把空间谱精度带来的影响减少到最小。取空间方向网格划分数N=256, 使得空间谱精度的误差可以忽略。为比较不同参数组合下的误差，选择两个参数A和ε的不同组合, 计算终止时间取t=1。为方便收敛阶测试, 时间步长从0.1开始依次减半, 相变量u的数值解的‖·‖2误差及收敛阶数的计算结果见表1。

表1 不同参数下时间方向的‖·‖2误差及阶数Table 1 The‖·‖2 errors and orders in temporal direction with different parameters

从表1结果可以看出, 尽管理论上对系数A的要求极高, 在数值模拟时却没有这一限制条件。无论是A=1,A=10还是A=20,都没有满足定理3中的条件, 但是计算结果的精度非常高。另一方面, 随着A值不断增加, 绝对误差越来越大, 这是因为A越大会产生更大的数值耗散。

对于空间方向的误差估计, 选择A=1,ε=0.025以及时间步长Δt=10-4。 时间步长取这么小的原因也是不希望时间层上的误差给空间误差的验证带来过多的干扰。 计算终止时间t=1，误差结果见表2。显然，Fourier谱精度随着空间网格划分数N的增加快速提升, 但是后期无论N值怎么增加, 精度都没有大幅提高, 因为此时已经不是空间离散误差而是时间离散误差占据主要作用。

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

3.2 晶体相场方程的动力学数值仿真

设二维空间区域为Ω=[0，128]×[0，128], 数值格式式(14)- 式(16)中的参数选取为A=2,ε=0.025, 空间步长固定为h=1, 晶体相的初始状态值为：

ui,j=0.06+0.03ri,j

(35)

其中ri,j为[0,1]服从均匀分布的随机数。为提高计算效率, 采取不同的时间步长进行数值模拟, 即在时间区间[0,100]上取Δt=0.01, 在[100,2 000]上取Δt=0.1, 晶体相场的动态演化在不同时间的状态见图1。

图1 不同时刻晶体相场的状态图Fig.1 State diagram of the crystal phase field at different times

图2给出不同时刻下的离散能量和离散质量的计算结果。 可以看出离散能量的变化和理论一致, 即保持单调不增。前一段由于时间步长较小, 使得图2中能量变化和时间步长较大的后一段相比较为平缓；中间段能量的变化较剧烈, 但随着能量衰减到一定程度, 能量变化呈现出平缓状, 此时晶体的相或者密度分布基本处于平稳状态。质量随着时间的变化，始终保持恒定，也验证了定理1的正确性。

图2 晶体离散能量和离散质量的动态演化图Fig.2 Dynamic evolution of the crystal discrete energy and discrete mass

4 结论

本文结合晶体相场模型的能量泛函结构, 利用凸分裂思想, 对周期边界条件的晶体相场方程提出了一种具有能量稳定性的数值算法。算法在空间上采用Fourier拟谱逼近，在时间采用上后向差分的三阶精度逼近, 并增加Duglas-Dupont正则项以保证数值格式的能量稳定性。对数值格式解的存在唯一性以及格式的能量稳定性进行了证明，采用数值仿真验证了格式在时间上的三阶精度和空间上谱精度, 并给出了方程的长时间动力学演化模拟。