陳慧琴,楊斌鑫

(太原科技大學(xué)應(yīng)用科學(xué)學(xué)院，太原030024)

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

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

1 相場(chǎng)方程的建立

1.1 相場(chǎng)控制方程

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

在相場(chǎng)方程演化過(guò)程中，假設(shè)了與時(shí)間有關(guān)的Ginzburg-Landau或Allen-Cahn方程：

(1)

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

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

(2)

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

其中局部自由能的表達(dá)式為：

flocal(φ,m)=

(3)

flocal(φ,m)對(duì)φ的一階導(dǎo)數(shù)為：

(4)

梯度能量的表達(dá)式為：

(5)

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

而ε的值可通過(guò)以下表達(dá)式計(jì)算：

(6)

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

(7)

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

(8)

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

(9)

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

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

(10)

1.2 溫度場(chǎng)控制方程

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

(11)

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

2 用傅里葉譜方法演化相場(chǎng)模型

2.1 傅里葉譜方法

傅里葉空間中空間導(dǎo)數(shù)的一般關(guān)系如下[11]：

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

特別的，當(dāng)n=1時(shí)且在二維空間中有[12]：

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

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

空間一階偏導(dǎo)：

時(shí)間一階偏導(dǎo)：

對(duì)于時(shí)間的偏導(dǎo)數(shù)采用向前差分。

拉普拉斯算子：

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

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

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

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

2.2 相場(chǎng)模型的離散

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

(12)

(12)的半隱式形式是：

(13)

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

(14)

對(duì)溫度場(chǎng)方程(11)進(jìn)行傅里葉譜變換：

(15)

(15)的半隱式形式是：

(16)

Δt是時(shí)間步長(zhǎng)n+1和n之間的時(shí)間增量，通過(guò)重新排列(16)，得：

(17)

3 結(jié)果與討論

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

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

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

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

圖1 初始晶核

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

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

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

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

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

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

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

圖5 加噪聲后枝晶凝固的時(shí)間演變

4 總結(jié)

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