王 崗，劉 艷，徐宗暢，肖 崗，張耀予，譚 凱，林一歆,

?

晶粒正常生長的Monte Carlo模擬

王 崗1，劉 艷1，徐宗暢1，肖 崗2，張耀予1，譚 凱2，林一歆1,2

1華中科技大學中歐清潔與可再生能源學院，武漢 4300742華中科技大學能源與動力工程學院，武漢 430074

本文建立了晶界能各向同性情況下晶粒生長的二維Monte Carlo模型，并對等溫情況下的晶粒生長過程進行了模擬。在模擬過程中，對傳統(tǒng)Monte Carlo方法中能量與概率統(tǒng)計方法進行了改進。為了更加直觀地顯示出晶粒生長過程中系統(tǒng)能量的變化，統(tǒng)計了在整個晶粒生長過程中能量的變化趨勢，結果與晶粒尺寸變化相符合。模擬得到的晶粒生長指數在0.35 ~ 0.45之間，與理論值相符，證明了改進方法的可靠性。

計算機模擬；晶粒生長；Monte Carlo方法

Monte Carlo (MC) 方法是研究晶粒生長以及其相關過程的一種有效方法。使用MC方法研究晶粒生長的模型實際上是從Ising模型和Potts模型演化出來的，這兩種方法最初被用于研究鐵磁體系統(tǒng)。使用MC方法研究晶粒生長實際上就是將不同取向的晶體賦予不同的取向數，根據概率性原則來決定晶粒生長這一傳播過程。

最早將MC方法運用到微觀結構演化過程的人是Anderson等人[1]。后來，Anderson等人又將這種方法應用到晶粒生長[2]、含有雜質粒子時的晶粒生長[3]，異常晶粒生長[4-6]等過程。隨后，Q-state Potts模型開始被用于研究兩相晶粒生長的機理[7,8]、晶粒生長過程中晶體結構的變化[9]以及外部條件對晶粒尺寸分布的影響[10]等。近年來，關于晶粒生長的MC模擬的研究主要著眼于對傳統(tǒng)晶粒生長模型的修改以使其更加貼合實際的晶粒生長過程，例如：Mason等人[11]對模型中Hmiltonian函數以及轉換概率進行了修改，模擬了動力學各向異性晶粒生長；Fang等人[12]運用新的模擬算法模擬了單相和兩相陶瓷材料在燒結過程中材料組織結構的演化；Wang等人[13]通過改進的MC算法模擬了晶體的初次再結晶過程和退火孿晶現象；Allen等人[14]運用統(tǒng)計學方法，量化晶粒各向異性的程度，然后通過MC方法模擬了單相材料的各向異性生長。

國內學者對晶粒生長MC模擬的研究也取得了一些顯著的進展。秦湘閣等人[15]運用MC方法對單相多晶體各向同性顆粒組織進行了三維建模并模擬其生長過程；鐘曉征等人[16,17]運用MC方法模擬了多晶材料晶粒生長，并模擬了正常晶粒生長和異常晶粒生長兩種情況；張繼祥等人[18]提出了一種“擇優(yōu)轉換”方法對MC方法進行了改進，模擬了晶粒的生長過程；王海東等人[19]根據晶粒生長機理建立改進的轉換概率模型模擬焙燒過程中的晶粒生長，并對不同溫度，不同焙燒時間以及不同激活能條件下的晶粒生長過程進行了模擬；王浩等人[20]運用MC方法對3種現存的三維個體晶粒生長速率拓撲依賴性方程進行了仿真驗證；馬非等人[21]采用改進的MC模擬方法對不同燒結溫度和保溫時間下的晶粒生長演化過程進行計算機模擬，并使用多孔Al2O3陶瓷燒結實驗過程中的晶粒生長數據對模型進行了驗證。

在實際應用中，傳統(tǒng)的MC方法存在效率低下的問題，這是因為在傳統(tǒng)MC方法中，格點取向改變范圍是在1 ~之間取一隨機數作為轉變取向。在實際的晶粒生長過程中，晶粒生長實際上只受其相鄰格點影響，其余格點對其取向轉變影響微乎其微。如果將所有格點取向都考慮進去，難免會造成計算效率低下的情況。因而，本文根據實際情況對晶粒取向轉變的算法進行了改進，在格點重新定向時，不再考慮所有的取向值，而僅考慮其鄰近格點的取向值，這樣就大大縮減了計算時間，提高了計算的準確率。

1生長動力學

1952年，Burke發(fā)現晶粒的平均尺寸與晶粒生長時間之間存在一定的關系：

(1)

式中，1是一個與溫度和晶界移動率有關的常數，0是系統(tǒng)在初始時刻的平均晶粒尺寸，是晶粒生長指數。關于晶粒生長指數的取值歷來都有爭議，因為在不同實驗條件、不同模型條件下模擬出的結果都不相同，一般而言，在不考慮二相粒子以及晶粒純度對于晶粒生長的影響的條件下，值一般趨近于2。

晶粒生長速率與晶粒尺寸之間的關系是由Hillert[22]提出的，其關系式如下：

(2)

其中，是一個與晶界遷移率和晶界能有關的常數，是單個晶粒的半徑，cr是臨界半徑。

晶界全都是高能量區(qū)域，晶界的存在會增加系統(tǒng)的自由能，從而使系統(tǒng)處于一種亞穩(wěn)定的狀態(tài)[23]。當外界提供一定的能量波動 (例如退火過程)，系統(tǒng)中的晶界就會通過自身調節(jié)，使晶界處于新的(亞) 穩(wěn)定狀態(tài)，從而降低整個系統(tǒng)的能量。

按具體形貌討論總的晶界能可以更簡單地討論每一個晶界對結合點的作用力，從而獲得晶界結合點的平衡條件。如圖1所示，晶界結合點的平衡條件可以表示為：

(3) 圖1晶界穩(wěn)定條件Figure 1 Stable grain boundary conditions

晶界張力在三叉點處的平衡使晶界發(fā)生彎曲。在彎曲界面上存在一個指向界面曲率中心的力，作用在單位面積上的力用表示，在這里稱為晶界移動驅動力，大小由Laplace's方程給出：

(4)

其中，1和2為晶界的兩個主曲率半徑，為界面張力。在二維系統(tǒng)中，2?￥，故上式可簡化為：

(5)

設在驅動力作用下，晶界移動一段距離，則能量變化D為

(6)

其中為晶界遷移方向與方向的夾角。

假設系統(tǒng)中界面張力為常數，由于驅動力的作用，晶粒長大過程中晶界將始終朝其曲率中心方向移動，所以具有內凹型晶界的晶粒將長大，晶界外凸的晶粒將縮小，而平直晶界因其曲率半徑無限大，驅動力等于零，將靜止不動。

圖2所示為一系列等邊長的邊型晶粒。由圖2可見，大于六邊的結構具有內凹邊界，他們在退火過程中會不斷長大；而小于六邊的為外凸晶界，在退火過程中不斷減小。這就是所謂的-6規(guī)則。Mullins[24]指出晶體面積的變化率與其拓撲關系存在如下的線性關系：

(7)

其中是指定晶粒的邊界數，是一個與晶界遷移率和晶界能成比例的常數。式 (7) 表明，當晶粒邊界數大于6時，晶粒長大；相反，當小于6時，晶粒將收縮；當邊界數等于6時，晶粒將不生長。

圖2晶粒拓撲結構

Figure 2 Grain topological structures

2模擬方法

MC-Potts模型是一種離散的統(tǒng)計方法。因為它十分適用于多相系統(tǒng) (例如固相與液相)，所以它被廣泛應用于擴散、離子傳輸、相變等現象的模擬研究。

利用Potts模型模擬晶粒生長時，首先應把晶粒形貌離散化，如圖3所示[25]。每一個離散格點都被賦予1 ~之間的數值，這個數值對應于晶粒在微觀下的取向。研究表明，的取值應當大于32[25]，因為只有這樣才能減少取向數對于模擬結果的影響。

Figure 3 Discrete structure

晶界的能量由不同格點之間的相互作用給出，用Hmiltonian方程來描述：

(8)

式中，是克羅內克符號，是相鄰格點之間的作用力，是周圍格點總數 (本文我們取如圖3所示點陣，因此任意一個格點周圍有8個鄰接格點)，q表示選定格點取向，q表示相鄰格點的取向。

在MC-Potts模型模擬晶粒生長的過程中，我們首先隨機選取一個格點，計算該格點取向與相鄰格點是否發(fā)生轉變，轉變概率由式 (9) 給出：

(9)

其中，D表示取向改變過程中能量的變化，由式 (10c) 計算得出，是玻爾茲曼常數，是絕對溫度。

在傳統(tǒng)的MC方法中，能量變化是這樣計算的：首先隨機選取一個格點，已知這個格點的取向值為，那么，我們就在除了以外的1 ~個取向中隨機選取一個取向，根據下列公式計算出該格點取向變?yōu)橄噜徣∠蚝?，系統(tǒng)能量差值的變化D，進而求出其轉換概率：

(10a) (10b) (10c)

但是在實際的晶粒生長過程中，我們所選取的格點生長方向實際上只與其相鄰晶粒有關。傳統(tǒng)MC方法將所有的晶粒取向作為格點可能將要轉變的方向顯然是不合適的。在本文中，對于所選取的格點，其轉變取向必然要在其相鄰8個格點中選取。為了更符合實際情況，我們將這8個格點的取向都加以考慮，分別計算出取向轉變后的能量改變 (D1~D8)，比較這8個能量差值的大小，選取絕對值最大的D(即取向轉變后系統(tǒng)能量損失最大) 對應的取向作為轉變取向，再根據式 (9) 求出其轉換概率。

3 模擬流程

對晶粒生長的Monte Carlo模擬流程如下 (圖4)：

(1) 輸入模擬初始數據 (取向數、格點數、模擬步數MCS等)；

(2) 建立晶粒初始形貌，將其離散到格點上；

(3) 對每一個格點隨機賦予取向數，代表該格點處的晶粒取向；

(4) 隨機選取一個格點，根據式 (8) 計算格點能量；

(5) 計算該格點取向轉變?yōu)橄噜徣∠蚝蟮哪芰浚?/p>

(6) 計算結點取向數改變前后能量差Δ；

(7) 根據式 (9) 計算轉換概率;

(8) 在模擬過程中產生[0,1] 之間的隨機數，比較和，判斷取向改變是否被接受：w3時接受；<時不接受。

(9) 如果取向改變被接受，則該格點取向數改變?yōu)槠湎噜徃顸c，晶粒長大；反之，如果改變未被接受，該晶粒則沒有長大；

(10) 一個MCS步完成；

(11) 再重新選取一個格點，按照(4) ~ (9) 的步驟模擬，直至完成所有循環(huán)。

按照上述模擬流程使用Matlab軟件編寫了軟件，并對模擬結果進行了后處理。模擬采用的點陣為100′100，取向值設定為32，模擬步數MCS為1000。

Figure 4 Flow diagram of the simulation process

4模擬結果與分析

圖5所示為模擬建立的晶粒初始形貌，圖6則給出了晶粒生長過程中晶粒形貌隨時間的演化情況?？梢钥闯?，在= 0 MCS ~ 250 MCS之間，晶粒生長迅速，尺寸變化明顯，大晶粒吞并小晶粒，晶界交點處多為120度，證明該模擬很好地符合了實驗現象。而在250 MCS后，晶粒尺寸變化并不明顯，于是我們統(tǒng)計了晶粒平均尺寸等參數來進行比較分析。

圖5晶粒初始形貌 (= 0)

Figure 5 The initial grain morphology (= 0)

晶粒平均尺寸的變化規(guī)律如圖7所示?？梢钥闯觯Я？焖偕L發(fā)生在0 MCS ~ 200 MCS這個階段；進入400 MCS后，晶粒尺寸基本保持穩(wěn)定。

圖7晶粒平均尺寸隨模擬時間的變化
Figure 7 Variation of average grain size with simulation time

為了進一步研究晶粒生長過程，我們統(tǒng)計了每一個MCS步時整個晶粒所儲存的能量，以便清楚地顯示出晶粒長大過程中整個系統(tǒng)的能量變化，結果如圖8所示。可以看出，晶粒生長過程實際上就是晶粒晶界儲能減少的過程，隨著晶粒儲能的減少，晶粒生長。在0 MCS ~ 200 MCS之間儲能急劇下降，對應于圖7中晶粒尺寸的急劇上升；400 MCS后，晶粒尺寸保持穩(wěn)定，晶粒儲能基本不變。

圖8 系統(tǒng)總能量隨模擬時間的改變

Figure 8 Variation of total system energy with simulation time

根據模擬結果，我們可以得到晶粒生長曲線如圖9所示：晶粒生長指數在0 MCS ~ 200 MCS之間有一個較大的變化，但最后穩(wěn)定在0.35 ~ 0.45之間，這與理論值0.5差別不大。研究表明，理論值0.5只是針對理想狀況，實際上由于雜質、間隙等因素的影響，大部分材料晶粒生長指數應當在0.1 ~ 0.4之間，這與我們的模擬結果相符。

圖9晶粒生長指數隨時間的變化

Figure 9 Variation of grain growth exponent with simulation time

5結論

本文根據傳統(tǒng)晶粒生長的Monte Carlo模型，提出了格點重新取向過程中新的隨機提取方法。模擬分析證明了該改進算法的可靠性。將晶粒平均尺寸與系統(tǒng)能量相比較，驗證了晶粒長大過程對應于系統(tǒng)能量減少這一現象。模擬得到的晶粒生長指數在0.35 ~ 0.45之間，與使用傳統(tǒng)Monte Carlo方法模擬得到的結果 (= 1/3[26]) 相比，更加貼近于理論結果。

[1] ANDERSON MP, SROLOVITZ DJ, GREST GS, et al. Computer simulation of grain growth: I, kinetics [J]. Acta Metallurgica, 1984, 32 (5): 783-791.

[2] SROLOVITZ DJ, ANDERSON MP, SAHNI PS, et al. Computer simulation of grain growth: II, grain size distribution, topology, and local dynamics [J]. Acta Metallurgica, 1984, 32 (5): 793-802.

[3] SROLOVITZ DJ, ANDERSON MP, GREST GS, et al. Computer simulation of grain growth: III, influence of a particle dispersion [J]. Acta Metallurgica, 1984, 32 (9): 1429-1438.

[4] SROLOVITZ DJ, GREST GS, ANDERSON MP. Computer simulation of grain growth: V, abnormal grain growth [J]. Acta Metallurgica, 1985, 33 (12): 2233-2247.

[5] ROLLETT AD, SROLOVITZ DJ, ANDERSON MP. Simulation and theory of abnormal grain growth: isotropic grain boundary energies and mobilities [J]. Acta Metallurgica, 1989, 37 (4): 1227-1240.

[6] ROLLETT AD, MULLINS WW. On the growth of abnormal grains [J]. Scripta mETALLURGICA ET mATERIALIA, 1996, 36 (9): 975-980.

[7] OKUDA K and rollett AD. Monte Carlo simulation of elongated recrystallized grains in steels [J]. Computational Materials Science, 2005, 34 (3): 264-273.

[8] COUTURIER G, DOHERTY R, MAURICE C, et al. 3D finite element simulation of the inhibition of normal grain growth by particles [J]. Acta Materialia, 2005, 53 (4): 977-989.

[9] IVASISHIN OM, SHEVCHENKO SV, VASILIEV NL, et al. 3D Monte-Carlo simulation of texture- controlled grain growth [J]. Acta Materialia, 2003, 51 (4): 1019-1034.

[10] ZOLLNER D, STREITENBERGER P, RIOS PR. Shedding some light on the early grain growth regime: about the effect of the initial microstructure on normal grain growth [J]. Computational Materials Science, 2016, 113: 11-20.

[11] MASON JK, LIND J, LI SF, et al. Kinetics and anisotropy of the Monte Carlo model of grain growth [J]. Acta Materialia, 2015, 82: 155-166.

[12] FANG B, HUANG C, XU C, et al. Numerical simulation of microstructural evolution of ceramic tool materials [J]. Computational Materials Science, 2008, 44 (2): 707-715.

[13] WANG W, HELBERY AL, BRISSET F, et al. Monte Carlo simulation of primary recrystallization and annealing twinning [J]. Acta Materialia, 2014, 81: 457-468.

[14] ALLEN JB, CORNWELL CF, DEVINE BD, et al. Simulations of anisotropic grain growth in single phase materials using Q-state Monte Carlo [J]. Computational Materials Science, 2013, 71: 25-32.

[15] 秦湘閣, 劉國權. 多晶體晶粒尺度三維組織建模及可視化[J]. 北京科技大學學報, 2001, 23 (6): 519-522.

[16] 鐘曉征, 陳偉元, 王豪才, 等. 多晶材料晶粒生長的Monte Carlo計算機模擬方法：Ⅰ, 模擬正常晶粒生長[J]. 功能材料, 1999, 17 (3): 9-12.

[17] 鐘曉征, 陳偉元, 王豪才, 等. 多晶材料晶粒生長的Monte Carlo計算機模擬方法：Ⅱ，模擬異常晶粒生長[J]. 功能材料, 999, 17 (3): 13-15.

[18] 張繼祥, 關小軍, 孫勝. 一種改進的晶粒長大Monte Carlo模擬方法[J]. 金屬學報, 2004, (5): 457-461.

[19] 王海東, 張海, 李海亮, 等. 焙燒過程晶粒生長的Monte Carlo模擬[J]. 中國有色金屬學報, 2007, 17 (6): 990-996.

[20] 王浩, 劉國權, 秦湘閣. 三維晶粒長大速率方程的大尺度Potts模型Monte Carlo仿真驗證[J]. 金屬學報, 2008, 44 (1): 13-18.

[21] 馬非, 申倩倩, 賈虎生, 等. 用改進的Monte Carlo算法模擬多孔Al2O-3陶瓷燒結過程中的晶粒生長[J].人工晶體學報, 2011, (5): 1299-1304.

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

[23] 毛衛(wèi)民, 趙新兵. 金屬的再結晶與晶粒長大[M]. 北京: 冶金工業(yè)出版社, 1994.

[24] MULLINS WW. Two-dimensional motion of idealized grain boundaries [J]. Journal of Applied Physics, 1956, 27 (8): 900-904.

[25] ZHANG D, WENG G, GONG S, et al. Computer simulation of grain growth of intermediate- and final-stage sintering and Ostwald ripening of BaTiO3-based PTCR ceramics [J]. Materials Science and Engineering B, 2003, 99 (1): 428-432.

[26] 葉日晴, 趙建華, 何陵輝. 退火過程中晶粒生長的二維計算機模擬[J]. 無機材料學報, 2001, (1): 122-128.

Monte Carlo Method for Normal Grain Growth Similation

WANG Gang1, LIU Yan1, XU Zong-Chang1, XIAO Gang2, ZHANG Yao-Yu1, TAN Kai2, LIN Yi-Xin1,2

1China-EU Institute for Clean and Renewable Energy, Huazhong University of Science＆Technology, Wuhan 430074, China2School of Energy and Power Engineering, Huazhong University of Science＆Technology, Wuhan 430074, China

A two-dimensional Monte Carlo model for the grain growth in under the isotropic condition was established and used to simulate the normal grain growth process. In the simulation process, some changes in the part of energy calculation and probability statistics in the traditional Monte Carlo method were made. The simulation results confirmed that the energy evolution of the system is corresponding to the grain size change. The obtained grain growth exponent,= 0.35 ~ 0.45, is close to its theoretic value, proving the reliability of the improved method.

Computer simulation; Grain growth; Monte Carlo method

TB34

1005-1198 (2016) 06-0434-08

A

10.16253/j.cnki.37-1226/tq.2016.07.002

2016-06-20

2016-10-31

國家自然科學基金 (21203069)。

王 崗 (1991-), 男, 河南湯陰人, 碩士研究生。E-mail: wanggangty44@gmail.com。

林一歆 (1982-), 女, 湖北武漢人, 副教授。E-mail: yixinlin@hust.edu.cn。