晶粒生長的雙區域動力學蒙特卡洛模擬

張　平,　李建華,　李洪林

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

?

晶粒生長的雙區域動力學蒙特卡洛模擬

張平,李建華,李洪林

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

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

動力學蒙特卡洛; 晶粒生長; 微觀組織; 晶粒生長指數

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

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

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

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

1　晶粒生長的kMC模擬

1.1晶粒生長

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

(1)

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

(2)

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

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

Δθ=

(3)

其中:SG1和SG2為晶界兩側格點內的取向數;Q為總取向數。

1.2晶粒生長的傳統動力學蒙特卡洛方法

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

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

(4)

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

(2)每個事件(網格點)有相應的發生轉換的概率ρ,所有事件組成一個概率事件列表,

(5)

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

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

(5) 返回步驟(2)。

1.3格點能量計算的傳統方法及其改進

傳統kMC方法中計算格點能量公式如下:

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

(6)

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

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

(7)

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

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

為了進一步模擬真實的晶粒形貌的演化,本文將格點分為晶界區和非晶界區,并將勢函數引入到晶界區格點能量計算。考慮到不同取向的晶粒之間存在不同能量,用該格點能量代替傳統的計算結果實現模擬。由于只在晶界區引入勢函數,所以減少了計算量。

2　晶粒生長的雙區域kMC模擬方法

2.1現有kMC計算晶界區格點能量的缺點

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

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

2.2模擬晶粒生長的雙區域kMC方法

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

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

圖2　區域劃分示意圖Fig.2　Schematic diagram of regional division

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

2.2.3非晶界區域的計算非晶界區域的生長方法如上述傳統kMC算法所示,對非晶界區域的格點隨機賦予取向數,再通過計算格點的轉換概率來判斷是否產生新的取向。

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

(8)

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

(9)

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

(10)

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

格點取向數的轉換概率計算公式如下：

(11)

最后再通過計算得到的取向數的轉換概率判斷得到格點取向。

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

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

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

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

3　實驗分析

3.1實驗設置

實驗模擬采用spparks開源軟件[12]實現了一個原型系統來模擬晶粒生長。參數設置為50×50網格,格點取向數為180種,勢函數采用LJ勢函數。

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

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

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

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

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

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

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

3.3雙區域kMC與傳統kMC的對比

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

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

圖7、圖8分別示出了傳統動力學蒙特卡洛和雙區域動力學蒙特卡洛的模擬結果。將兩種kMC方法對比發現,采用雙區域kMC的模擬方法得到的模擬結果更加接近理論值。

圖7　傳統動力學蒙特卡洛模擬結果Fig.7　Simulation results of traditional kinetic Monte Carlo

圖8　雙區域動力學蒙特卡洛模擬結果Fig.8　Simulation results of double-region kinetic Monte Carlo

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

3.4雙區域kMC與改進kMC的對比

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

圖9　近鄰轉換機制Fig.9　Nearest neighbor conversion mechanism

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

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

圖10　晶粒生長指數Fig.10　Grain growth exponent

4　結束語

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

[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]張繼祥,關小軍.異常晶粒長大的Monte Carlo模擬[J].中國有色金屬學報,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]張繼祥,關小軍,孫勝，等.晶粒長大過程微觀組織演變Monte Carlo方法模擬[J].山東大學學報(工學版),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-),男,江蘇常州人,碩士生,研究方向為計算機輔助設計。E-mail:zhangping_1@foxmail.com

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

TP391