晶粒生長的雙區(qū)域動力學(xué)蒙特卡洛模擬

張　平,　李建華,　李洪林

(華東理工大學(xué)信息科學(xué)與工程學(xué)院,上海 200237)

?

晶粒生長的雙區(qū)域動力學(xué)蒙特卡洛模擬

張平,李建華,李洪林

(華東理工大學(xué)信息科學(xué)與工程學(xué)院,上海 200237)

現(xiàn)有的動力學(xué)蒙特卡洛方法模擬晶粒生長過程的結(jié)果與理論計算值之間大多存在偏差,其中重要的原因是現(xiàn)有方法大多采用只適用于計算小角度晶界Read-Shockley公式。針對此問題,本文提出了一種對大、小角度晶界都適用的雙區(qū)域動力學(xué)蒙特卡洛模擬方法。該方法通過基于勢函數(shù)的格點能量獲得取向概率,從而避免使用Read-Shockley公式。為了解決引入勢函數(shù)導(dǎo)致計算量大的問題,劃分動力學(xué)蒙特卡洛的單區(qū)域結(jié)構(gòu)為晶界區(qū)域和非晶界區(qū)域的雙區(qū)域結(jié)構(gòu),并將勢函數(shù)的計算限制在晶界區(qū)內(nèi),從而減少了計算時間。實驗?zāi)M結(jié)果表明:與傳統(tǒng)的和改進的動力學(xué)蒙特卡洛方法相比,雙區(qū)域動力學(xué)蒙特卡洛方法的微觀組織演變更符合晶粒生長特征,晶粒生長指數(shù)也更接近理論值。

動力學(xué)蒙特卡洛; 晶粒生長; 微觀組織; 晶粒生長指數(shù)

晶粒生長是影響多晶材料物理和機械性能的一個重要過程。目前,已有很多種理論和實驗方法用于研究晶粒生長,其中比較常見的方法有相場法[1]、蒙特卡洛(Monte Carlo)法[2]和動力學(xué)蒙特卡洛(Kinetic Monte Carlo,kMC)法[3-6]。

動力學(xué)蒙特卡洛方法是用于模擬自然界中某些隨時間演變過程的一種蒙特卡洛方法[7]。由于動力學(xué)蒙特卡洛方法的基礎(chǔ)是隨機性,使得它非常適合描述某些過程,例如空位擴散、生長、化學(xué)反應(yīng)等。

眾多學(xué)者在用動力學(xué)蒙特卡洛方法模擬晶粒生長時均使用Read-Shockley 模型[8]計算晶界能。該模型由Read和Shockley于1950年提出[9],由于其只適用于小角度,對大角度晶界的計算結(jié)果往往偏離真實值。國內(nèi)外的學(xué)者對此都提出了改進方法。張繼祥等[10]提出將大角度晶界的能量設(shè)為常數(shù);Mallick等[11]在計算小角度晶界能時用Read-Shockley模型,在計算大角度晶界能時嘗試將公式進行標準化擴展,雖然擴展公式的變化趨勢符合原公式,但是計算得到的晶界區(qū)的能量卻與原來的值存在差別。雖然已有學(xué)者對Read-Shockley模型提出了改進,但是改進方法對于大角度晶界的計算仍然無法得到理想值。

為此,本文提出將勢函數(shù)引入晶界區(qū)格點能量的計算中,避免Read-Shockley模型的局限性,模擬了各向異性條件下晶粒生長的組織結(jié)構(gòu)演化過程。

1　晶粒生長的kMC模擬

1.1晶粒生長

多晶材料中的晶粒生長是一個自發(fā)的過程,可以形象地形容為“大晶粒吞并小晶粒”。由于實際晶粒組織的復(fù)雜性,對晶粒的形狀簡化假設(shè)都是基于統(tǒng)計平均的方法。Hillert[12]從統(tǒng)計角度出發(fā),提出晶粒長大模型如下:

(1)

其中:M為晶界遷移率;γ為晶界區(qū)的能量;R為晶粒半徑;Rc為臨界半徑;α為常數(shù)。實驗研究[13]和計算機模擬[14]結(jié)果表明,晶界區(qū)的能量可以改變晶粒的生長行為。傳統(tǒng)的計算晶界區(qū)能量的方法是采用Read-Shockley公式,見式(2)。

(2)

其中:θ′為晶粒間取向差,常設(shè)置為0到π的某個值;θ為晶界夾角;J為大角度晶界能。

公式中的實際晶粒取向通過如下方法與模擬中隨機生成的取向數(shù)相對應(yīng)[10]:在直角坐標系中,對于二維中的晶粒有0～360°的取向,由于晶粒取向的特殊性可表示為0～180°,因此用180個取向數(shù)表示,每個取向數(shù)代表一個對應(yīng)的晶體取向。晶界兩側(cè)的取向差Δθ(即晶界夾角)的計算如下:

Δθ=

(3)

其中:SG1和SG2為晶界兩側(cè)格點內(nèi)的取向數(shù);Q為總?cè)∠驍?shù)。

1.2晶粒生長的傳統(tǒng)動力學(xué)蒙特卡洛方法

在介觀尺度觀察晶粒生長的過程中,可以將該過程看成是晶粒從一種狀態(tài)到另一種狀態(tài)的躍遷,而kMC是適用于系統(tǒng)狀態(tài)轉(zhuǎn)變的方法。kMC是由Bortz等[15]在蒙特卡洛算法基礎(chǔ)上提出的NFW(N-Fold Way)方法演化而來,目前已有十幾種不同版本的kMC方法。

kMC方法的具體思路是:首先確定系統(tǒng)中發(fā)生的事件的個數(shù)M以及每個事件發(fā)生變化的速率kij,然后由[0,1]之間的隨機數(shù)確定從初始態(tài)到下一個狀態(tài)發(fā)生變化所需的時間τ,以及選擇要發(fā)生的事件,發(fā)生相應(yīng)的變化后系統(tǒng)將從初始態(tài)躍遷到下一個狀態(tài),如此反復(fù),直到所設(shè)定的限制條件。躍遷時間的計算公式如下:

(4)

在晶粒生長中,動力學(xué)蒙特卡洛方法首先將多晶體材料的連續(xù)微觀組織結(jié)構(gòu)映射到離散的二維網(wǎng)格上,即對于每個網(wǎng)格點隨機賦予一個整數(shù)值,具有相同整數(shù)值的相鄰網(wǎng)格點表示一個晶粒,不同整數(shù)的交界作為晶界。具體的生長算法如下:

(2)每個事件(網(wǎng)格點)有相應(yīng)的發(fā)生轉(zhuǎn)換的概率ρ,所有事件組成一個概率事件列表,

(5)

其中:常數(shù)γmax和Mmax分別表示晶界移動的最大能量和遷移率;ΔE表示晶界區(qū)格點能量;Ts為模擬溫度。

(3)依次選取每個事件嘗試新取向,計算概率,然后在[0,1]內(nèi)隨機產(chǎn)生一個分布均勻的數(shù)ε，ε=ερ0，如果ρ>ε，則表示該網(wǎng)格點的新取向是可以被接受的,也就表示重新取向可被接受,則晶界發(fā)生遷移;如果ρ≤ε,則表示該網(wǎng)格點的新取向是不可以被接受的,也就是此網(wǎng)格點會繼續(xù)保持原來的取向。

(5) 返回步驟(2)。

1.3格點能量計算的傳統(tǒng)方法及其改進

傳統(tǒng)kMC方法中計算格點能量公式如下:

本系統(tǒng)設(shè)計了左右兩個測試子系統(tǒng)進行交替測試。一個測試，一個準備，以保證整個測試過程中測試時間的0等待。測試流程包括以下四個測試子過程：1.裝置在放置區(qū)和測試區(qū)之間進出的自動控制，2.是裝置硬件接點的自動檢測；3.是針對非數(shù)字化裝置進行的交流頭精度自動校準及標簽的自動打?。?.左右測試子系統(tǒng)的自動切換。整機測試流程的智能控制意在解決以上四個過程的自動控制問題，見圖4。

(6)

其中:γ(θij)表示在格點間隨著取向差θij變化的單位面積的能量，γ(θij)計算采用Read-Shockley公式，即公式(2);δSiSj表示Kronecker函數(shù);n為近鄰的網(wǎng)格點數(shù)。

目前已有的改進kMC中比較典型的方法是Mallick[11]提出的將Read-Shockly中大角度計算進行標準化擴展,擴展公式如下:

(7)

其中:r為常數(shù),取0.618;Δθ為取向差。圖1示出了采用傳統(tǒng)和擴展的kMC計算得到的大角度晶界能。傳統(tǒng)kMC中大角度晶界能設(shè)置為常數(shù),而擴展的kMC大角度晶界能隨取向差的增大而逐漸增大。但是擴展的kMC中,當(dāng)取向差角度大于20°時,晶界能成為固定值,這與實際情況不符合。

圖1　擴展Read-Shockley公式擬合大角度晶界能Fig.1　Extended Read-Shockley formula fitting energy in large angle grain boundary

為了進一步模擬真實的晶粒形貌的演化,本文將格點分為晶界區(qū)和非晶界區(qū),并將勢函數(shù)引入到晶界區(qū)格點能量計算?？紤]到不同取向的晶粒之間存在不同能量,用該格點能量代替?zhèn)鹘y(tǒng)的計算結(jié)果實現(xiàn)模擬。由于只在晶界區(qū)引入勢函數(shù),所以減少了計算量。

2　晶粒生長的雙區(qū)域kMC模擬方法

2.1現(xiàn)有kMC計算晶界區(qū)格點能量的缺點

在各向同性模擬中,晶界區(qū)格點能量可視為相同的值;在各向異性模擬中,則必須將其設(shè)為不同值,即將晶界分為大角度和小角度兩種晶界。而當(dāng)Read-Shockley公式用于各向異性的計算時,由于其只適用于小角度晶界的計算,對于大角度的計算結(jié)果會出現(xiàn)偏差[16],因此會導(dǎo)致最終的模擬生長過程偏離預(yù)測值。本文在原有方法中將勢函數(shù)加入晶界區(qū)格點能量計算中,重新計算得到取向數(shù)的轉(zhuǎn)換概率,避免使用Read-Shockley公式。

晶粒的長大依賴于晶界區(qū)提供的驅(qū)動力。引入勢函數(shù)計算時,需要對每個格點分別進行計算,增加了計算量,因此本文只在晶界區(qū)引入勢函數(shù)進行計算。

2.2模擬晶粒生長的雙區(qū)域kMC方法

2.2.1概述雙區(qū)域kMC方法主要分為非晶界區(qū)和晶界區(qū),分別計算各自的取向數(shù)。在晶界區(qū)域中引入勢函數(shù),通過計算格點能量引入到取向數(shù)轉(zhuǎn)換概率的計算中,優(yōu)點是不用區(qū)別大、小角度晶界,而且可以對不同類型的晶界使用合適的勢函數(shù),以此來彌補傳統(tǒng)方法的不足。由于kMC格點取向數(shù)的改變會影響周圍取向數(shù)的變化,因此兩區(qū)域kMC能夠較好地連接。

2.2.2非晶界區(qū)域和晶界區(qū)域的劃分如圖2所示,整個區(qū)域的劃分依據(jù)是勢函數(shù)的截斷半徑。由于采用的是Lennard-Jones勢函數(shù),所以將距離交界線的兩個格點的長度作為晶界區(qū)域,其余為非晶界區(qū)域。圖中三角形表示晶界區(qū)域,橢圓表示非晶界區(qū)域,不同的顏色表示界線兩側(cè)不同取向數(shù)的晶粒。

圖2　區(qū)域劃分示意圖Fig.2　Schematic diagram of regional division

劃分晶界區(qū)域的大小對生長情況和計算量存在著影響。相同模擬時間內(nèi),如果晶界區(qū)擴大,則計算量增大,但是晶粒平均尺寸相同,原因是超過截斷半徑的晶粒之間無法相互影響;相同模擬時間內(nèi),如果晶界區(qū)縮小,則計算量減小,但是晶粒平均尺寸減小,原因是驅(qū)動晶粒生長的晶界能減小。

2.2.3非晶界區(qū)域的計算非晶界區(qū)域的生長方法如上述傳統(tǒng)kMC算法所示,對非晶界區(qū)域的格點隨機賦予取向數(shù),再通過計算格點的轉(zhuǎn)換概率來判斷是否產(chǎn)生新的取向。

2.2.4晶界區(qū)域的計算多晶材料是由許多取向不同的單晶組成的,不同種取向之間的區(qū)域稱為晶界區(qū)。由于同種取向間的作用力小于不同種取向間作用力,因此在晶界區(qū)就有取向不同而產(chǎn)生的力,稱為取向力。采用勢函數(shù)直接計算取向力,然后再求晶界區(qū)格點能量。勢函數(shù)選用Lennard-Jones勢函數(shù),LJ勢函數(shù)方法簡單、計算量小,公式如下:

(8)

其中:ε表示兩個原子間相互吸引作用的強弱;σ為作用勢等于0時原子間的距離。

(9)

模擬溫度的計算公式如下：

(10)

其中:s為取向力;v=s×d距離,d距離表示力方向的移動距離。

格點取向數(shù)的轉(zhuǎn)換概率計算公式如下：

(11)

最后再通過計算得到的取向數(shù)的轉(zhuǎn)換概率判斷得到格點取向。

2.2.5非晶界區(qū)域和晶界區(qū)域的連接kMC方法中的取向數(shù)是向二維網(wǎng)格中八邊形的點陣賦予一個整數(shù)值(如圖3所示)。

在系統(tǒng)的取向數(shù)變化過程中,當(dāng)相鄰格點間取向一致時,則式(9)中的δSiSj為0,不一致時則為1。因此kMC方法中的每個格點取向數(shù)的計算會影響到其周圍鄰居點的取向數(shù)的變化,即晶界區(qū)域取向數(shù)的改變會導(dǎo)致非晶界區(qū)域的變化。

圖3　模擬體系中四邊形點陣Fig.3　Quadrilateral lattice in simulation system

在雙區(qū)域kMC方法中,通過只在晶界區(qū)中引入勢函數(shù)計算得到晶界區(qū)格點能量,然后通過取向數(shù)的計算將兩區(qū)域較好地連接,避免了由Read-Shockley公式引起的不足,模擬結(jié)果也更符合理論規(guī)則。

3　實驗分析

3.1實驗設(shè)置

實驗?zāi)M采用spparks開源軟件[12]實現(xiàn)了一個原型系統(tǒng)來模擬晶粒生長。參數(shù)設(shè)置為50×50網(wǎng)格,格點取向數(shù)為180種,勢函數(shù)采用LJ勢函數(shù)。

3.2各向異性對晶粒長大特征的影響

3.2.1組織演變圖4所示為晶粒的組織結(jié)構(gòu)隨模擬時間的演化圖,采用了各向異性條件。圖中不同類型的晶粒用不同的顏色表示。可以看出,在各向異性條件下,組織中晶粒尺寸趨于非均勻,形狀趨于不規(guī)則晶粒,說明各向異性阻礙晶粒長大,降低組織均勻化程度。

圖4　各向異性條件下雙區(qū)域kMC模擬的晶粒 生長微觀組織演變Fig.4　Microstructures evolution in anisotropic grain growth with bi-region kMC

3.2.2晶粒尺寸分布平均晶粒尺寸的測量是在模擬二維圖中隨機布置3條平行于邊界的直線(橫向或縱向),計算每條直線上的晶粒平均截距(覆蓋在一個晶粒上的測量線段的長度稱為截距),取算術(shù)平均值作為晶粒平均直徑[17]。其中,晶粒平均截距=長度/截到的晶粒個數(shù)。

本文隨機選取縱向的3條直線,統(tǒng)計直線所截的晶粒尺寸和平均尺寸,并計算算術(shù)平均值。

圖5示出了最大晶粒直徑與晶粒平均直徑的比值Db/Dm,MCS表示蒙特卡洛時間步長。圖中顯示,Db/Dm比值大于2.5,但小于4。以上模擬結(jié)果符合各向異性條件下的正常晶粒長大特征,與文獻[18]相符。

圖5　最大晶粒尺寸與平均尺寸比值Fig.5　Maximum grain size and average size ratio

3.3雙區(qū)域kMC與傳統(tǒng)kMC的對比

圖6示出了不同晶界夾角情況下,由Read-Shockley公式計算得到的晶界區(qū)能量變化過程?？梢钥闯?雖然變化趨勢基本符合要求,但是對于大角度晶界(即夾角大于15°)不符合實際情況[19]。從圖中也發(fā)現(xiàn),生長過程中大角度晶界所占的數(shù)量多于小角度晶界。

圖6　Read-Shockley公式擬合晶界區(qū)格點能量Fig.6　Read-Shockley formula fitting energy in grain boundary

圖7、圖8分別示出了傳統(tǒng)動力學(xué)蒙特卡洛和雙區(qū)域動力學(xué)蒙特卡洛的模擬結(jié)果。將兩種kMC方法對比發(fā)現(xiàn),采用雙區(qū)域kMC的模擬方法得到的模擬結(jié)果更加接近理論值。

圖7　傳統(tǒng)動力學(xué)蒙特卡洛模擬結(jié)果Fig.7　Simulation results of traditional kinetic Monte Carlo

圖8　雙區(qū)域動力學(xué)蒙特卡洛模擬結(jié)果Fig.8　Simulation results of double-region kinetic Monte Carlo

晶粒生長過程中主要有兩種拓撲演化機制:小于六邊形晶粒的消失和近鄰切換機制。計算晶粒的邊數(shù)可以通過計算其周圍不同晶粒的種類獲得。

3.4雙區(qū)域kMC與改進kMC的對比

本文選取改進kMC中比較典型的擴展Read-Shockly方法[11]與雙區(qū)域kMC方法做實驗比較。通過圖9和圖10可以發(fā)現(xiàn)雙區(qū)域方法比改進kMC更接近理論值。

圖9　近鄰轉(zhuǎn)換機制Fig.9　Nearest neighbor conversion mechanism

圖9展示的是用雙區(qū)域kMC模擬的近鄰轉(zhuǎn)換機制,而用改進的kMC無法得到該結(jié)果。近鄰轉(zhuǎn)換機制指的是由于在晶粒生長過程中,晶界能量會不斷減小,因此相鄰的兩個晶粒交接處的晶界會消失。如圖9中,4和5這兩個相鄰的晶粒之間的交界線消失。

圖10示出了晶粒生長指數(shù)比較結(jié)果。從圖中可以看出，改進的kMC的晶粒生長指數(shù)超過了0.51,超過了理論最大值(晶粒長大指數(shù)n與Hillert[15]指數(shù)的關(guān)系n=1/m,根據(jù)統(tǒng)計方法,理論最大值n=0.5),而雙區(qū)域kMC符合各向異性的晶粒生長指數(shù)的變化規(guī)律,即小于各向同性條件下的生長速度,也遠小于最大理論值。

圖10　晶粒生長指數(shù)Fig.10　Grain growth exponent

4　結(jié)束語

本文介紹了利用雙區(qū)域動力學(xué)蒙特卡洛方法對各向異性的晶粒生長進行模擬,解決了傳統(tǒng)蒙特卡洛模擬只適用于計算小角度晶界的限制。該方法將勢函數(shù)通過計算格點能量引入取向數(shù)的轉(zhuǎn)換概率計算中,并采用分區(qū)域方法減少計算量。雙區(qū)域動力學(xué)蒙特卡洛方法與已有的傳統(tǒng)動力學(xué)蒙特卡洛方法和改進的動力學(xué)蒙特卡洛方法的模擬結(jié)果相比,在相同的模擬時間內(nèi),雙區(qū)域方法的微觀組織演變更符合晶粒生長特征,晶粒生長指數(shù)也更接近理論值。

[1]FAN D,CHEN L Q.Computer simulation of grain growth using a continuum field model[J].Acta Materialia,1997,45(2):611-622.

[2]YU Q,NOSONOVSKY M,ESCHE S K.On the accuracy of Monte Carlo Potts models for grain growth[J].Journal of Computational Methods in Sciences & Engineering,2008,8(4/6):227-243.

[3]FU K,FU Y.Kinetic Monte Carlo study of metal organic chemical vapor deposition growth mechanism of GaSb quantum dots[J].Applied Physics Letters,2008,93(10):101906-10190-3.

[4]ANDERSON M P,SROLOVITZ D J,GREST G S,etal.Computer simulation of grain growth：I.Kinetics[J].Acta Met allurgica,1984,32(5):783-791.

[5]PATTERSON B R,ROWENHORST D J,TIKARE V,etal.Affinities for topological arrangements in grain structures[J].Acta Materialia,2014,79(30):411-420.

[6]HOLM E A,HASSOLD G N,MIODOWNIK M A.On misorientation distribution evolution during anisotropic grain growth[J].Acta Materialia,2001,49(15):2981-2991.

[7]ROHRER G S.Influence of interface anisotropy on grain growth and coarsening[J].Annual Review of Materials Research,2005,35(1):

[8]YU Q,ESCHE S K.Modeling of grain growth kinetics with Read-Shockley grain boundary energy by a modified Monte Carlo algorithm[J].Materials Letters,2002,56:47-52.

[9]BATTAILE C C.The kinetic Monte Carlo method:Foundation,implementation,and application[J].Computer Methods in Applied Mechanics & Engineering,2008,197:3386-3398.

[10]張繼祥,關(guān)小軍.異常晶粒長大的Monte Carlo模擬[J].中國有色金屬學(xué)報,2006,16(10):1689-1697.

[11]MALLICK A,VEDANTAM S.Phase field study of the effect of grain boundary energy anisotropy on grain growth[J].Computational Materials Science,2009,46(1):21-25.

[12]HILLERT M.On the theory of normal and abnormal grain growth[J].Acta Met Allurgica,1965,13(3):227-238.

[13]GOTTSTEIN G,SHVINDLERMAN L S,GOTTSTEIN G,etal.Grain boundary migration in metals:Thermodynamics,kinetics,applications[J].Journal of Engineering Mechanics,1998,126:888.

[14]GREST G S,SROLOVITZ D J,ANDERSON M P.Computer simulation of grain growth：Ⅳ.Anisotropic grain boundary energies[J].Acta Met allurgica,1985,33(85):509-520.

[15]BORTZ A B,KALOS M H,LEBOWITZ J L.A new algorithm for Monte Carlo simulation of king spin systems[J].Journal of Computational Physics,1975,17:10-18.

[16]GRUBER J,MILLER H M,HOFFMANN T D,etal.Misorientation texture development during grain growth：Part Ⅰ.Simulation and experiment[J].Acta Materialia,2009,57(20):6102-6112.

[17]RADHAKRISHNAN B,ZACHARIA T.Simulation of curvature-driven grain growth by using a modified monte carlo algorithm[J].Met allurgical & Materials Transactions A,1995,26(1):167-180.

[18]張繼祥,關(guān)小軍,孫勝，等.晶粒長大過程微觀組織演變Monte Carlo方法模擬[J].山東大學(xué)學(xué)報(工學(xué)版),2005,35(4):1-5.

[19]OLMSTED D L,FOILES S M,HOLM E A.Survey of computed grain boundary properties in face-centered cubic metals:I.Grain boundary energy[J].Acta Materialia,2009,57(13):3694-3703.

[20]MULLINS W W.Two-dimensional motion of idealized grain boundaries[J].J Appl Phys,1956,27:900-904.

Bi-region Kinetic Monte Carlo Simulation of Grain Growth

ZHANG Ping,LI Jian-hua,LI Hong-lin

(School of Information Science and Engineering,East China University of Science and Technology,Shanghai 200237,China)

There exist deviations between the results of most kinetic Monte Carlo (kMC) simulation of grain growth and the theoretical value of grain growth.An important reason is that the Read-Shockley formulation in most kMC simulation is only suitable for the computation of small-angle grain boundary.Aiming at this problem,a bi-region kMC method is proposed in this paper,which handles both situation of small-angle grain boundary and large-angle grain boundary.In the proposed method,the transition probability of orientation number is obtained using lattice energy based on potential functions,not using Read-Shockley formula.In order to solve the large computational problem caused by introducing potential functions,kMC single-region structure of traditional dynamics is divided into two regions:the one without grain boundary and the one with the grain boundary.Due to potential function computation only used in the region of grain boundary,the computing time is reduced.Experiment results show that the micro-structure evolution of bi-region kMC method is more consistent with the feature of grain growth than that of existing traditional and improved kMC,and the grain growth exponent is closer to the theoretical value than that of existing kMC.

kMC; grain growth; micro-structure; grain growth exponent

A

1006-3080(2016)03-0387-06

10.14135/j.cnki.1006-3080.2016.03.015

2015-09-25

張平(1991-),男,江蘇常州人,碩士生,研究方向為計算機輔助設(shè)計。E-mail:zhangping_1@foxmail.com

通信聯(lián)系人：李建華,E-mail:jhli@ecust.edu.cn

TP391