一類相場方程的能量穩定性分析及數值模擬

霍俊蓉, 劉 昊, 溫學兵, 張榮培, 蔚喜軍

(1. 沈陽師范大學 數學與系統科學學院, 沈陽 110034； 2. 廣東工業大學 應用數學學院, 廣州 510006； 3. 北京應用物理與計算數學研究所, 北京 100088)

0 引 言

隨著計算機技術的發展以及材料熱力學和動力學數據的完備, 相場模型逐漸成為一種模擬各種材料微結構演化過程的重要工具. 該數值模擬方法在生物和材料力學等領域應用廣泛, 如薄膜和裂紋擴展等. 在相場模型中, 需引入另一種連續的場變量以表示不同的相, 稱為相場. 任意成分或結構不同的多相或多域的微觀結構可通過一系列連續的場變量u表示,u隨時間和空間的演化過程可由偏微分方程描述. Cahn-Hilliard方程[1]是常見的相場模型方程, 該方程在物理、 化學、 材料科學[2-3]、 熱力學[4-5]、 流體力學[6-7]和圖像處理[8-9]等領域應用廣泛.

在不同物質相分離的過程中, 二元混合體系的2個組分自發分離, 并在各組分中形成純疇域. 根據Ginzburg-Landau理論, 二元混合體系的自由能方程為

(1)

(2)

(3)

該邊界條件表示區域邊界上的能量流為零.

由于Cahn-Hilliard方程為四階非線性擴散方程, 不易求得解析解, 因此, 一般借助數值方法求其近似解. 人們利用有限元方法[10-11]、 間斷有限元方法[12-15]、 Lattice Boltzmann方法[16]、 最小二乘譜元法[17]和譜方法[18-19]對Cahn-Hilliard方程進行了研究. 此外, Fourier多網格算法[20]、 局部加密純無網格有限點集法[21]和快速顯式算子分裂法[22]也被應用于求解該類方程. 該類方程空間離散后得到一組剛性非線性常微分方程組, 顯式時間離散方法對時間步長的限制較嚴格, 隱式時間離散方法允許較大的時間步長, 但需求解大型的非線性代數方程組. 本文應用有限差分方法空間求解具有恒定遷移率的Cahn-Hilliard方程, 利用Kronecker積寫出二維Laplace算子的微分矩陣并將其對角化. 結合快速離散余弦變換(fast discrete cosine transform, FDCT)實現快速求解. 該方法可證明全離散格式下的離散能量保持能量穩定性, 具有存儲量小和計算速度快等優點.

1 數值方法

在熱力學中, 相場模型通過數值求解演化方程, 得到場變量在模擬區域和時間節點上的值, 進而反映材料體系的微結構演化過程[23]. 場變量u的演化方程通常為非線性偏微分方程, 本文采用有限差分法以及Crank-Nicolson方法, 以變量離散取值后對應的函數值近似微分方程中獨立變量的連續取值, 進而求解Cahn-Hilliard方程, 該方法在解決相場模型問題的同時可保證解的穩定性.

1.1 有限差分方法

考慮遷移率M(u)=1的Cahn-Hilliard方程, 則式(2)等價于

(4)

網格步長為hx=(a-b)/Nx,hy=(c-d)/Ny.在每個網格節點處有其特定的場變量.設u在網格節點(xi,yj)的數值解為ui,j(1≤i≤Nx, 1≤j≤Ny), 定義離散解ui,j為Nx×Ny階的矩陣U.由齊次Neumann邊界條件, 以及網格內部點(xi,yj)和邊界處的二階導數差分格式, 得到該方程在x,y方向的微分矩陣分別為

將解矩陣按列向量化后得

U=vec(U)=(u1,1…uNx,1,u1,2…uNx,2,…,u1,Ny…uNx,Ny)T.

(5)

設Ix,Iy分別為Nx,Ny階單位矩陣, 利用微分矩陣Bx和By以及Kronecker積的定義[24], 可將式(4)的中心差分格式寫為

(6)

引理1對分別為Nx×Nx階與Ny×Ny階的正交矩陣Cx與Cy, 矩陣Bx與By有如下的特征分解:

(7)

(8)

引理2對N階矩陣X以及N階矩陣A,B,C,D, Kronecker積有如下性質[26]:

1) (A?B)(C?D)=AC?BD;

2) (A?B)-1=A-1?B-1;

3) (BT?A)vec(X)=vec(AXB);

4) (A?B)T=AT?BT.

(9)

引理3K為對稱負定矩陣, 且可分解為K=-LTL.

由引理3可得:

引理4對?U∈Nx×Ny, 有(KU,U)=-‖LU‖2, 其中‖·‖為Euclid范數.

1.2 Crank-Nicolson方法

Cahn-Hilliard方程為四階非線性方程, 顯式時間離散方法對時間步長的限制較嚴格, 隱式格式具有較好的穩定性. Crank-Nicolson方法是常用的二階隱式差分格式, 一般用于求解具有守恒型的非線性方程. 下面應用Grank-Nicolson差分方法對空間離散后得到的式(6)進行時間離散, 得

(10)

式(10)可等價為

(11)

其中

由式(1), 定義離散自由能形式為

(12)

將式(10)中的兩式分別與Tn+1/2和Un+1-Un做內積, 可得

(Un+1-Un,Tn+1/2)=(KTn+1/2,Tn+1/2)=-‖LTn+1/2‖2,

(13)

(Tn+1/2,Un+1-Un)=(-γKUn+1/2,Un+1-Un)+(Fn+1/2,Un+1-Un).

(14)

對Tn+1/2=-γKUn+1/2+Fn+1/2兩端與Un+1-Un做內積, 可得

(Tn+1/2,Un+1-Un)=(-γKUn+1/2,Un+1-Un)+(Fn+1/2,Un+1-Un).

(15)

對式(15)右端兩式分別進行運算, 有

則可將式(15)化為

(18)

由(Un+1-Un,Tn+1/2)=-‖LTn+1/2‖2以及二元混合體系的自由能ε(u)的定義, 可得

εn+1(U)-εn(U)≤0.

(19)

由式(19)可知, 在零通量邊界條件下, 有εn+1(U)≤εn(U).即離散形式的自由能具有關于時間非增的性質.

由于相場模型滿足一種被稱為能量穩定的非線性穩定關系, 即自由能泛函隨時間遞減.因此, 本文的數值格式可在快速求解的同時滿足離散能量穩定.

1.3 快速實現

下面應用FDCT快速求解全離散后的Cahn-Hilliard方程, 即式(11). 利用數值方法結合初始值與邊界條件, 求解與微分方程對應的離散數值解, 進而描述二元合金的凝固、 固態相變、 電子遷移、 粗化與晶粒生長等物理現象的演化過程.

首先對其進行整理得

(20)

令C=Cy?Cx, 可將矩陣K和K2對角化, 得到K=-CTΛC,K2=CTΛ2C, 結合式(14), 可將式(20)寫為

(21)

分三步求解式(21).

1) 計算CUn與CFn+1/2, 由引理2的性質3)可得

(22)

式(22)可由MATLAB中命令dct2實現.這里向量G和H分別由矩陣G和H向量化得到, 即G=vec(G),H=vec(H).

(23)

實現.

3) 計算CTG2和CTH2, 利用引理2中的性質2)和性質3), 可得

(24)

式(24)可由MATLAB中命令idct2實現.非線性代數方程組(20)可寫為

Un+1=G3-H3.

(25)

當

(26)

時, 迭代結束.

2 數值實驗

求解具有初始條件u0(x,y)=0.1sin(x)sin(y)的Cahn-Hilliard方程(4), 并驗證其自由能關于時間是非遞增的.

將其網格剖分為64×64, 時間離散步長為0.01, 過渡區域厚度γ=0.01.離散形式質量與能量隨時間的變化曲線如圖1所示.由圖1可見, 在全離散格式下, 質量守恒且自由能關于時間是非遞增的.因此, 本文采用的數值方法可在離散后保持Cahn-Hilliard方程的本質特征, 并保證全離散格式解的穩定性.

圖1 離散形式質量(A)與能量(B)隨時間的變化曲線Fig.1 Change curves of discrete mass (A) and energy (B) with time

為更好展示二元合金中金屬部件濃度的相分離過程及其粗化過程, 在[0,2π]2求解區域中求解初始條件為

的算例, 其中u為二元合金中金屬部件濃度.將求解區域網格剖分為128×128, 取計算時間T=2, 時間步長Δt=10-4, 金屬部件濃度過渡區域的參數γ=0.01.不同時間下的二元合金中金屬部件濃度如圖2所示.由圖2可見, 濃度變化包含2個階段.第一個階段為相位分離, 通常由自由能εn(U)的最小化引起, 使濃度在0

圖2 不同時間下的二元合金中金屬部件濃度Fig.2 Concentration of metal parts in binary alloys at different time

為驗證本文的數值方法可在保證穩定性的同時提高計算效率, 將網格剖分成不同大小, 對是否采用FDCT的計算時間進行比較. 取時間步長Δt=10-4, 最終計算時間T=0.005.在不同網格剖分條件下, 采用FDCT和未采用FDCT的CPU時間列于表1. 由表1可見, 本文結合FDCT使計算效率得到極大提高.

表1 采用FDCT和未采用FDCT的CPU時間比較

綜上, 本文在空間和時間上分別應用有限差分方法和Crank-Nicolson方法對具有恒定遷移率的Cahn-Hilliard方程進行全離散. 利用Kronecker積寫出二維Laplace算子的微分矩陣, 并引入Kronecker積的性質將微分矩陣進行特征分解. 在實現過程中, 采用不動點迭代法結合快速離散余弦變換求解, 以提高計算效率. 本文結合微分矩陣的性質證明了全離散格式下的Cahn-Hilliard方程仍保持其離散能量不增, 通過數值實驗驗證了該方法的穩定性.