一類相場方程的能量穩(wěn)定性分析及數(shù)值模擬

霍俊蓉, 劉 昊, 溫學(xué)兵, 張榮培, 蔚喜軍

(1. 沈陽師范大學(xué) 數(shù)學(xué)與系統(tǒng)科學(xué)學(xué)院, 沈陽 110034； 2. 廣東工業(yè)大學(xué) 應(yīng)用數(shù)學(xué)學(xué)院, 廣州 510006； 3. 北京應(yīng)用物理與計(jì)算數(shù)學(xué)研究所, 北京 100088)

0 引 言

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

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

(1)

(2)

(3)

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

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

1 數(shù)值方法

在熱力學(xué)中, 相場模型通過數(shù)值求解演化方程, 得到場變量在模擬區(qū)域和時(shí)間節(jié)點(diǎn)上的值, 進(jìn)而反映材料體系的微結(jié)構(gòu)演化過程[23]. 場變量u的演化方程通常為非線性偏微分方程, 本文采用有限差分法以及Crank-Nicolson方法, 以變量離散取值后對(duì)應(yīng)的函數(shù)值近似微分方程中獨(dú)立變量的連續(xù)取值, 進(jìn)而求解Cahn-Hilliard方程, 該方法在解決相場模型問題的同時(shí)可保證解的穩(wěn)定性.

1.1 有限差分方法

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

(4)

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

將解矩陣按列向量化后得

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

(5)

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

(6)

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

(7)

(8)

引理2對(duì)N階矩陣X以及N階矩陣A,B,C,D, Kronecker積有如下性質(zhì)[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為對(duì)稱負(fù)定矩陣, 且可分解為K=-LTL.

由引理3可得:

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

1.2 Crank-Nicolson方法

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

(10)

式(10)可等價(jià)為

(11)

其中

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

(12)

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

(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)

對(duì)Tn+1/2=-γKUn+1/2+Fn+1/2兩端與Un+1-Un做內(nèi)積, 可得

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

(15)

對(duì)式(15)右端兩式分別進(jìn)行運(yùn)算, 有

則可將式(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).即離散形式的自由能具有關(guān)于時(shí)間非增的性質(zhì).

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

1.3 快速實(shí)現(xiàn)

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

首先對(duì)其進(jìn)行整理得

(20)

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

(21)

分三步求解式(21).

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

(22)

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

(23)

實(shí)現(xiàn).

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

(24)

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

Un+1=G3-H3.

(25)

當(dāng)

(26)

時(shí), 迭代結(jié)束.

2 數(shù)值實(shí)驗(yàn)

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

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

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

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

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

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

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

表1 采用FDCT和未采用FDCT的CPU時(shí)間比較

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