枝晶生長的傅里葉譜相場模擬

陳慧琴,楊斌鑫

(太原科技大學應用科學學院，太原030024)

幾個世紀以來，液相凝固過程中復雜微觀結(jié)構的形成，如雪花的形成和金屬鑄態(tài)微觀結(jié)構，一直以來都是科學家們所關注的。尤其枝晶在凝固過程中的微觀結(jié)構尺度的演化決定了金屬的許多物理力學性質(zhì)[1]，是金屬實際加工中具有精確尺寸和優(yōu)良性能的重要依據(jù)和基礎科學問題，因為幾乎所有的金屬體系都起源于液態(tài)。作為一種免于追蹤固-液界面，準確模擬液相凝固過程中枝晶界面形態(tài)的重要工具，相場方法(PF)深受國內(nèi)外廣大學者的追捧和廣泛使用。相場法是以Ginzburg-Landau理論為物理基礎，由van der Waals[2]所提出的擴散界面模型擴展而來，主要優(yōu)點在于，相場方法的運用使得計算枝晶生長過程中免于追蹤復雜的固-液兩相界面，通過使用實值0和1，分別代表液相和固相，用0～1之間表示固液界面，簡化了運算過程。Wheeler等人建立了合金的第一個相場模型[3-4]，稱為WBM模型。Kim等人[5-6]介紹了另一種模型，采用了薄界面極限，被稱為KKS模型。

在求解相場方程數(shù)值解的過程中，有限差分法是使用最為普遍的求解離散方法。然而，有限差分法的使用會離散出計算極為復雜的九點差分格式，增加了計算運行的時間。傅里葉譜方法作為一種求解偏微分方程的數(shù)值計算方法，將偏微分方程的解近似地展開成光滑函數(shù)的有限級數(shù)展開式，再由所此有限級數(shù)展開式與偏微分方程，求解得到有限級數(shù)展開式的系數(shù)方程組。使用傅里葉譜方法求解相場方程，不需要繼續(xù)計算復雜的九點差分格式的算子，并且對于與周期性邊界條件，傅里葉級數(shù)的計算更為便捷，譜方法的精度更為準確，相比有限差分方法，傅里葉譜方法的使用大大加快了求解相場方程的運算速度，提高了數(shù)值計算與數(shù)值模擬的效率。

1 相場方程的建立

1.1 相場控制方程

下面給出的模型基于Kobayashi[7]最著名的工作，其作為樹枝狀凝固的最早相位場模型之一。該模型包括兩個變量：一個是非守恒相場參數(shù)φ(r,t)，它在固體中取值1，在液體中取值0.另一個是溫度場T(r,t)，也隨著凝固的進展而進化。其中，r是空間位置，t是時間。

在相場方程演化過程中，假設了與時間有關的Ginzburg-Landau或Allen-Cahn方程：

(1)

其中τ是松弛時間，Χ(φ,m)是系統(tǒng)自由能函數(shù)。

之后，考慮Ginzburg-Landau型自由能[8]：

(2)

上述積分項中的fgrad(φ)是梯度能量，flocal(φ,m)是局部自由能。

其中局部自由能的表達式為：

flocal(φ,m)=

(3)

flocal(φ,m)對φ的一階導數(shù)為：

(4)

梯度能量的表達式為：

(5)

其中ε是各向異性梯度能量系數(shù)，它決定界面層厚度。為了引入界面各向異性，假定在界面處ε依賴于外法向量的方向[9]。

而ε的值可通過以下表達式計算：

(6)

σ(θ)=1+δcos(j(θ-θ0))

(7)

上式中δ是各向異性強度；j是各向異性模數(shù)；θ0是初始偏移角，取作常數(shù)。這里角度θ被定義為：

(8)

(4)式中的m，假定它取決于過冷度和溫度[7]，它的表達式為：

(9)

其中α是一個正常數(shù)，Teq是平衡溫度。

將(2)，(3)，(4)，(5)代入到(1)中，整理可得：

(10)

1.2 溫度場控制方程

由焓守恒定律導出的溫度場T(r,t)的演化方程：

(11)

上式中T無量綱化的溫度，所以特征冷卻溫度是零，平衡溫度是1[10].κ是一種無量綱潛熱，它與潛熱成正比，與冷卻強度成反比。為了簡單起見，無論在固相中還是液相中都將擴散常數(shù)設置為相同的。

2 用傅里葉譜方法演化相場模型

2.1 傅里葉譜方法

傅里葉空間中空間導數(shù)的一般關系如下[11]：

其中{·}k是括號內(nèi)量的傅里葉變換，k是k個傅立葉模式的系數(shù)。

特別的，當n=1時且在二維空間中有[12]：

其中F(A)=fft(A)，F(xiàn)-1(A)=ifft(A)，fft傅里葉正變換，ifft是傅里葉逆變換。

根據(jù)上述表達，相場模型中的空間一階偏導、時間一階偏導、拉普拉斯算子的傅里葉變換如下：

空間一階偏導：

時間一階偏導：

對于時間的偏導數(shù)采用向前差分。

拉普拉斯算子：

F(2φ)=-k2{φ}k，F(xiàn)(2T)=-k2{T}k

若使用有限差分九點格式離散拉普拉斯算子則有：

?2ψ(i,j,n)=ψxx(i,j,n)+ψyy(i,j,n)=

很直觀的就可以觀察到，傅里葉譜方法在運算上比有限差分簡單的多。

2.2 相場模型的離散

用傅里葉譜方法變換方程(10)兩邊[13]：

(12)

(12)的半隱式形式是：

(13)

Δt是時間步長n+1和n之間的時間增量[14]，通過重新排列(13)，得：

(14)

對溫度場方程(11)進行傅里葉譜變換：

(15)

(15)的半隱式形式是：

(16)

Δt是時間步長n+1和n之間的時間增量，通過重新排列(16)，得：

(17)

3 結(jié)果與討論

為了證明傅里葉譜方法可用于求解此類方程，對求解結(jié)果進行了數(shù)值模擬。

下面給出處理后方程所需要的相場模型參數(shù)的取值：

這些模擬是針對一個方形模擬單元進行的，該模擬單元Nx=Ny=300，網(wǎng)格步長為dx=dy=0.03，時間步長為Δt=0.000 1，時間迭代步數(shù)為4 000.

初始結(jié)晶核的大小是x2+y2<10.0，如圖1所示。

圖1 初始晶核

實驗所需的其他無量綱化參數(shù)在表1中已給出，利用Matlab進行編程，并用Paraview可視化得到下列圖像(如圖2所示)，分別是在時間迭代步數(shù)為400、1 400、2 400、3 400時的圖像，其中圖2第一行是枝晶凝固的時間演變，圖2第二行是溫度場的時間演變。在無噪聲的情況下，采用粗糙網(wǎng)格進行數(shù)值模擬，得到了定性的結(jié)果，與楊斌鑫等人的結(jié)果非常一致。

表1 模擬所需的無量綱化參數(shù)

張晨輝提到隨各向異性強度δ增加，枝晶尖端不再分叉，且在主枝上出現(xiàn)了側(cè)枝，主枝變細，優(yōu)先生長；而平衡溫度Teq增大,枝晶側(cè)枝的生長速度變快，側(cè)枝的尖端半徑在減小。通過改變模擬時的各向異性強度、平衡溫度等參數(shù)，得到了如圖3、圖4所示的實驗結(jié)果，其中a1-a3(時間迭代步數(shù)分別為900、1 900、2 900)以及b1-b3(時間迭代步數(shù)分別為600、1 600、2 600)是數(shù)值結(jié)果，而a4和b4是實驗結(jié)果，與實驗結(jié)果對比，可以看到數(shù)值結(jié)果與實驗結(jié)果非常一致。

圖2 各向異性模數(shù)j=6，各向異性強度ε=0.02時，枝晶凝固的時間演變和溫度場的演變

圖3 各向異性模數(shù)j=6，各向異性強度ε=0.06時，所獲得的片晶形態(tài)

圖4 平衡溫度Teq=0.9時枝晶凝固的時間演變和實驗結(jié)果

在圖5中(分別是在時間迭代步數(shù)為400、1 400、2 400、3 400時的圖像)，在界面上加了小噪聲，這樣更接近現(xiàn)實。采用了隨機數(shù)字序列arφ(1-φ)表示噪聲，其中a為噪聲幅值，此處取a=1，x均勻分布在[-0.5,0.5]這個區(qū)間上。對比圖2第一行發(fā)現(xiàn)，加了噪聲后的圖像主枝稍微粗壯了些，側(cè)枝數(shù)量減少，枝晶生長速度減緩。

圖5 加噪聲后枝晶凝固的時間演變

4 總結(jié)

通過建立傅里葉譜方法計算一般相場方程的演化求解步驟，使用傅里葉譜方法求解模擬枝晶生長的相場模型，簡化了方程中拉普拉斯算子的演化，提高了整體運算速度與實驗模擬效率，通過改變模擬參數(shù)，定性研究了晶體生長形態(tài)的變化，證明了使用傅里葉譜方法求解枝晶生長的相場法的有效性，是一種更為便捷的計算方法。