修正的晶体相场方程的无条件能量稳定数值格式

梁译泓 贾宏恩,2,*

(1.太原理工大学数学学院,晋中,030600;2.智能优化计算与区块链技术山西省重点实验室,晋中,030619)

1 引言

晶体相场方程(PFC)是由Elder 等引入相场变量模拟材料在原子空间尺度和扩散时间尺度上的微观结构时提出的[1,2].这种方程用途广泛,可以模拟晶粒增长、外延增长、裂纹扩展等行为.PFC 方程经Stefanovic 等加入弹性相互作用后发展为修正的晶体相场方程(MPFC)[3,4].该方程可以自洽地模拟大尺度上的快速弹性松弛.但文献[5]指出,MPFC 方程的原始自由能在某些时间区间上会随时间增长,因此需要引入一个伪能量并且证明它是耗散的.与PFC 方程相比,MPFC 方程引入了高阶时间导数,这给数值格式的构造带来了困难.对此,文献[6,7]提出了非线性且能量稳定的数值格式,文献[8]对MPFC 方程给出了LDG 方法,文献[9]提出了一个全离散的有限元格式,文献[10] 通过IEQ 方法提出了3 个线性的、无条件能量稳定的数值格式,文献[11,12]通过SAV 方法提出了线性的、无条件能量稳定的格式.

为了构造MPFC 方程的二阶线性格式,本文基于文献[13]中的拉格朗日乘子法,对空间六阶的非线性方程进行降阶处理.本文结构如下:第二节介绍MPFC 方程并对其能量进行修正;第三节应用拉格朗日乘子法构建一个二阶的数值格式,然后证明这个格式是质量守恒、唯一可解及无条件能量稳定的;第四节对这个二阶格式进行严格的误差估计;最后一节通过数值算例对理论分析进行验证.

2 MPFC 模型

Swift-Hohenberg 型自由能泛函为[1]

其中Ω⊂Rd(d=2,3),变量ϕ(x,t)表示密度场,常数∊∈(0,1).修正的相场晶体模型可以写成一个伪梯度流的形式:

其中常数α>0,β >0,迁移率常数M >0,化学势µ定义为

其中变量满足周期边界条件.方程(2.2)的初始条件为ϕ(x,0)=ϕ0,ϕt(x,0)=ψ0.

对式(2.2)两端积分得到常微分方程

引入一个新变量ψ(x,t):=ϕt(x,t),则.用算子(-∆)-1作用于式(2.2)两边后代入式(2.3),再与ϕt作内积,可得MPFC 方程满足下述能量耗散律:

3 二阶数值格式

通过引入两个辅助变量u(x,t)=ϕ2,w=∆ϕ,式(2.1)可以改写为

伪能量式(2.4)可以改写为

式(2.2),(2.3)的等价形式为

因为式(3.5)是一个常微分方程,故u不需要边界条件,其余变量取周期边界条件,初始条件为

用(-∆)-1作用于式(3.2)两边,再与ϕt作内积,并将式(3.6)代入得

然后将式(3.5)两边与ϕt,分别做内积得

再将式(3.6)两边与-w做内积得

最后,综合式(3.7)–(3.10),得

注3.1 对式(3.5)关于时间求积分可得u=ϕ2,把它和式(3.6)代入式(3.2)与(3.3)可得到式(2.2)和(2.3),因此,(2.2)与(3.2)–(3.6)这两个系统是等价的.(3.2)–(3.6)关于修正的能量式(3.1)是无条件稳定的.本文将对系统(3.2)–(3.6)构造数值格式,使其满足新能量律(3.11).

以下,用δt表示时间步长,tn=nδt(0≤n≤N=),ϕn表示ϕ(tn)的数值近似,C泛指某个正常数,它与δt和n无关,在不同情况下可以取不同的值.

下面对系统(3.2)–(3.6)构造二阶格式.

对于1≤n≤N-1,若已知(ϕn,ψn,un,wn) 和ϕn-1,则通过下述方式求解(ϕn+1,ψn+1,un+1,wn+1):

其中初始条件为ϕ0=ϕ0,ψ0=0,u0=(ϕ0)2.

注3.2 式(3.12)–(3.16)和式(3.17)–(3.21)都是通过引入拉格朗日乘子u=ϕ2得到的,它可以看作IEQ 方法[2]的一个特例.与文献[2]不同的是本文显式处理2∆ϕ,得到一个强制的双线性形式a(·,·).

对式(3.12)–(3.16)进行等价变形得到

其中,

由式(3.22)可得到ϕn+1,µn+1,wn+1,然后通过式(3.14)和式(3.15)可得到un+1,ψn+1.这样,引入的辅助变量u和ψ没有增加额外的计算成本.

由于对于具有周期边界条件的函数f和g都有

且(P(f),f)=0 当且仅当f≡0 时成立,故线性算子P(ϕ)是对称正定的.

对式(3.12)和(3.17)两边关于x积分得

对式(3.14)与(3.19)两边关于x积分得到

这说明系统满足质量守恒性.

下面讨论系统(3.22)(或其等价系统(3.12)–(3.16))的适定性.

定理3.1 线性系统(3.22)存在唯一解.

用(-∆)-1算子作用于式(3.23)两边,然后将式(3.24),(3.25)代入,可得

令V=,则式(3.26)的变分问题为:求ϕ∈V,使得

其中,

由（∇(-∆)-1ϕ,∇(-∆)-1ϕ）≤∥(-∆)-1ϕ∥·∥ϕ∥≤∥∇(-∆)-1ϕ∥·∥ϕ∥及Poincaré 不等式可得∥∇(-∆)-1ϕ∥≤C∥ϕ∥.因此,a(ϕ,v)=γ（(-∆)-1ϕ,v）+（P(ϕ),v）+≤C∥ϕ∥2∥v∥2,即a(·,·)是连续的.

由Lax-Milgram 定理知,变分问题(3.27)的解ϕ∈V存在且唯一,因而等价系统(3.22)的解存在且唯一.

类似地可证明初始格式(3.17)–(3.21)的解存在且唯一.定理证毕.

定理3.2 格式(3.12)–(3.16)是无条件能量稳定的,即

其中,

证明用算子(-∆)-1作用于式(3.12)两边,然后两边与ϕn+1-ϕn作内积,再结合式(3.14)可得

用ϕn+1-ϕn与式(3.13)两边作内积得

再将式(3.16)与wn=∆ϕn相减,然后两边与wn+1+wn作内积得

最后,结合式(3.29)–(3.32)可得

定理证毕.

定理3.3 格式(3.17)–(3.21)是无条件能量稳定的,即

其中,

证明用(-∆)-1算子作用于式(3.17)两边,然后与(ϕn+1-ϕn)作内积,再结合式(3.19)可得

将(ϕ1-ϕ0)与式(3.18)两边作内积,得

将式(3.21)与式w0=∆ϕ0相减,然后两边与wn+1做内积,得

最后,结合式(3.35)–(3.38)可得

定理证毕.

4 误差估计

本节对半离散格式(3.12)–(3.16)进行误差分析.为简单起见,用f≲g表示f≤Cg.根据离散的能量耗散律(3.28)和(3.34),对数值解ϕn有如下的H2估计.

定理4.1 假设初始条件足够正则,且≤C,则||ϕn||2≤C,其中0≤n≤N.

证明由式(3.34)得≤C.由式(3.39)得≤C.因此,

由式(3.33)得

故

所以,

又因为

所以,

由式(3.32)和式(3.38)知

进一步,有

定理得证.

设(3.2)–(3.6)的精确解(ϕ,ψ,u,w)满足下列正则性条件:

定义误差函数:

定理4.2 设(3.2)–(3.6)的精确解满足正则性条件(4.3).则有

证明将(3.2)–(3.6)改写为

其中,

它们满足

由式(4.5)与式(3.20)–(3.22),可得

用算子(-∆)-1作用于式(4.6)两边,然后与作内积,可得

结合式(4.11)–(4.13)得

上式右端各项的估计如下:

利用这些估计可得

定理证毕.

接下来对二阶格式(3.12)–(3.16)进行误差估计.

将式(3.2)–(3.6)改写为

其中,

对于截断误差式(4.15),有下述定理,其证明略.

定理4.3 在正则性条件式(4.3)下,截断误差式(4.15)满足

由式(4.14)及式(3.12)–(3.16),可得下列误差方程:

定理4.4 设式(3.2)–(3.6)的精确解(ϕ,ψ,u,w)满足正则性条件(4.3).则存在某个常数C1使得

证明用算子(-∆)-1作用于式(4.17)两边,然后与作内积得

下面分析式(4.23)右端各项.

结合式(4.22)–(4.28)得:

下面对式(4.29)的右端各项进行估计.

类似可得:

综上可得:

对式(4.30)两边令n取值1 到m,然后求和(1≤m≤N-1),可得

再由式(4.16)得

再由Gronwall 不等式知,存在某个常数C1,使得

定理证毕.

定理4.5 设(3.2)–(3.6)的精确解满足正则性条件式(4.3).则

即离散能量式(3.29)是对伪能量式(2.4)的二阶近似.

证明根据定理4.4 及三角不等式可得:

进一步,我们有

5 数值实验

本节的数值实验都在二维空间中进行,变量均使用周期边界条件,并且选择P1有限元空间.例5.1 计算数值格式(3.15)–(3.19)的时间精度,例5.2 模拟相变过程.

例5.1 方程参数设为α=1,β=0.9,M=1,∊=0.025,总时间T=1.0,积分区域Ω=[0,1]2,相变量的精确解为

为了与方程(2.2)保持一致,需要在方程(2.2)的右端增加一个外力项f(x,y,t),即

这里需要将ϕ(x,y,t)的精确值代入f(x,y,t)的表达式中进行计算.

变量ϕ,u的L2误差表达式分别为

变量ϕ,u关于时间的收敛阶的表达式分别为

实验结果如表1 所示.从表中可以看到随着时间步长减小,误差值越来越小,趋于0,收敛阶o(ϕ)与o(u)逐渐趋于2.由此说明本文所构造的差分格式在时间上是二阶收敛的,与文中理论分析的结论一致.

表1 二阶格式(3.15)–(3.19)的实验结果

例5.2 将积分区域Ω=(0,6π)×(0,6π)分成50×50 个正方形,每个正方形沿左下角(右上角) 对角线划分,得到一致三角剖分.方程参数设为α=0.1,β=1,M=1,∊=2,时间步长δt=0.25,初始条件为

图1 到图4 展示了液滴从条形分布到三角形分布的大致变化过程,实验进行到t=1250 时液滴的三角形分布达到稳定.

图1 t=5

图2 t=26.25

图3 t=30

图4 t=1250