基于物质跃迁能垒改进正常晶粒长大的Monte Carlo模拟

刘建元，张晓泳，周科朝

(中南大学 粉末冶金国家重点实验室，长沙 410083)

基于物质跃迁能垒改进正常晶粒长大的Monte Carlo模拟

刘建元，张晓泳，周科朝

(中南大学 粉末冶金国家重点实验室，长沙 410083)

在采用改进的Monte Carlo方法模拟正常晶粒的长大过程时，考虑到实际晶粒长大过程中物质状态跃迁需克服一定能垒，把跃迁能垒引入到状态跃迁概率计算中，使模拟过程所含物理意义更为明确。模拟结果表明：在各种模拟条件下得到的模拟晶粒长大指数为0.472～0.493，接近理论的晶粒长大指数(0.5)。改变模拟温度和能垒等模拟条件后，得到的晶粒形貌演变过程、晶界分布拓扑特征以及晶粒长大速率等模拟结果也均与晶粒长大动力学相关理论和实际晶粒长大规律相吻合。

Monte Carlo方法；晶粒长大模拟；晶粒长大指数；跃迁概率；物质跃迁能垒

Abstract:Considering that an energy barrier should be conquered to realize the mass transition during the real grain growth process, the Monte Carlo method was modified by introducing the energy barrier for mass status transition into the evaluation of mass fransition probability to simulate the normal grain growth. The aim of modification was to improve the physical significance of simulation. The simulation results show that the simulated grain growth exponent is 0.472−0.493, which is close to the theoretical value of 0.5. With the variation of simulation temperature and energy barrier, the simulation results about the grain morphology development, the topological characteristic of grain boundary and the variation of grain growth rate are all in accordance with the theories and real rules of normal grain growth kinetics.

Key words:Monte Carlo method; simulation of grain growth; grain growth exponent; transition probability; energy barrier for mass transition

Monte Carlo方法(简称MC)因具有模拟实现方法和程序设计简单、模拟速度相对较快以及能够利用计算机图形学直观描述微观结构等诸多优点，被广泛用于晶粒长大过程的模拟研究[1−6]。在MC模拟过程中，首先，建立离散化模拟对象，并赋予每个格点某一取向度值，由取向度值相同的相邻格点组成晶粒，由每对取向度不同的相邻格点组成晶界；然后，随机选取晶界处的格点，并尝试赋予其新取向度值，通过计算由此引起的晶界能变化，得到赋予被选格点新取向度值的概率，以此决定是否赋予被选格点新取向度值，进而实现物质在晶界处跃迁以及晶界迁移的模拟。对所有格点遍历处理一次即为一个模拟时间步(MCS)，并最终模拟出整个晶粒的长大过程。

在目前大多数文献所报道的MC模拟晶粒长大研究中，在给定模拟温度时，所得模拟结果与实际晶粒长大过程中的微观结构拓扑演变关系较为相符，且模拟的晶粒长大指数也接近理论长大指数(0.5)[6−9]。但是在模拟不同温度下的晶粒长大过程时，会出现模拟温度升高引起晶粒长大速率减慢这一与实际现象不符的模拟结果[10]。为此，文献[11−13]建立了如下的模拟时间—真实时间的线性和非线性对应关系：

式中：tMCS为模拟时间；T和t分别为真实温度和时间；A、K1和n1都为建模系数；K0为与原子振动频率相关的常数；Q为激活能；Z为单位晶界面积上原子的平均个数；Vm为原子体积；γ为晶界能；L0为初始平均晶粒尺寸；ΔSf为融化熵；λ为阵点距离；Na为阿伏伽德罗常数；R为摩尔气体常数；h为普朗克常数。利用上述对应关系可模拟出不同温度和升温速率下的正常晶粒长大过程。但在其模拟处理方式中，只涉及格点重新赋值后需降低晶界能这一要求，并没有考虑格点重新赋值要克服一定能垒以及格点通过能量起伏并克服能垒后有可能向高能态跃迁这两点，与实际晶粒长大过程中的物质状态跃迁在一定程度上不相符[14]。与此同时，针对原有MC模拟方法中跃迁概率计算存在的一些问题，文献[15−17]提出了新的跃迁概率计算表达式。但这些拟合公式缺乏较为真实的物理含义。如文献[16]提出了如下拟合公式：

式中：ΔE为格点重新赋值前、后晶界能的变化；TR为材料再结晶温度；a为拟合系数；K为波尔兹曼常数。

针对这些研究现状，本文作者首先分析原有的MC模拟方法中因模拟处理方式与实际物质跃迁过程不符而导致出现的一些问题，然后提出相应的改进措施，将物质在跃迁前需克服一定能垒这一因素在跃迁概率计算过程中直接体现，使模拟过程的物理含义更为明确，并用于模拟不同温度下的晶粒长大过程。

1 模拟过程

图1给出了MC模拟晶粒长大过程中用于描述微观结构格点集合的示意图。其中如果某一格点周围存在取向度值与之不同的其它相邻格点，则该格点位于晶界处，有可能被改变取向度值状态，进而引起晶粒长大和晶界迁移。

图1 MC模拟晶粒长大过程中用于描述微观结构格点集合的示意图Fig.1 Schematic of lattice integration representing microstructure in MC method to simulate grain growth

目前，大多数MC模拟晶粒长大研究中的跃迁模拟示意图如图2所示。对于图2中取向度值为2的中间被选格点i，首先赋予其周围与之不同的格点取向度值为1，并计算取向度值改变后引起的晶界能变化ΔE：

式中：J为晶界能；Qi和Qi′分别为被选格点i重新赋值前、后的取向度值；n为与格点i相邻的格点数；Qj为相邻格点的取向度值；δ为Kronecker算符(Qi或Qi′ ≠Qj时δ= 0；Qi或Qi′ =Qj时δ= 1)。若 ΔE≤0，则赋予格点i新取向度值Qj，即引起晶界处的物质迁移，并使总晶界能减小；若ΔE＞0，则计算格点i是否接受新取向度值Qj，进而得到引发高能态跃迁的Boltzman概率(P)，然后，在(0, 1)产生一个随机数W；当P≥W时，赋予格点i新取向度值Qi′，否则，仍保持原有取向度值。

式中：k为Boltzman常数；T为模拟温度。

通过分析可以发现，在上述跃迁模拟处理方式中，随着模拟温度T的升高，式(2)中物质向高能态跃迁的Boltzman概率相应增加。即在相同模拟时间下，随着kT的增大，所选格点被赋予新取向度值后引起晶界能增加(ΔE＞0)，晶界增多的模拟现象更容易发生，而发生格点改变取向度值后引起晶界能降低(ΔE≤0)这一事件的概率始终为 1，最终导致出现晶粒长大速率随模拟温度的升高而减慢这一与实际晶粒长大规律不符的模拟结果。在实际晶粒长大过程中，晶界处的物质须先克服一定能垒(G)后才能跃迁至新状态，进而实现晶界迁移。但这一因素在原有晶粒长大MC模拟方法中并未体现。为此，本研究提出如下改进：

图2 大多数MC模拟晶粒长大过程中物质状态跃迁的模拟示意图Fig.2 Schematic of mass status transition simulation used in previous MC method to simulate grain growth

1) 将模拟对象离散成n×n的四方格点，每个格点随机赋予某一整数型取向度值 Q(1≤Q≤Qmax)，由此完成模拟对象初始化，本研究中n=200，Qmax=160。

2) 随机选取任一处于晶界处的格点 i(即被选格点 i周围有取向度不同的相邻格点)，其取向度值为Qi，然后将每个相邻格点的取向度值Qj(Qi≠ Qj)随机逐一尝试赋予被选格点 i，并根据式(4)计算相应的晶界能变化ΔE，进而给出引起晶界能降幅最大(ΔE≤0)的取向度值Qjj。若将每个相邻格点的取向度值赋予被选格点i后，均引起晶界能增大(ΔE＞0)，则选取引起晶界能升幅最小的取向度值Qjj；若有多个取向度值赋予被选格点i后均引起晶界能降幅最大或者升幅最小，则从中随机选取。这种“择优转换”的处理方式考虑了体系能量降低的最大程度是体系状态优先变化的动力，并促使晶界朝其曲率中心移动，同时也可有效提高模拟效率[6]。

3) 根据下式给出是否赋予被选格点i新取向度值Qjj的概率(P)，然后在(0, 1)产生一个随机数W，当P≥W时，赋予格点i新取向度值Qjj，否则，仍保持原有取向度值。式中：G为通过改变被选格点取向度值来模拟物质跃迁时所需克服的激活能垒。

在式(6)给出的跃迁概率计算方法中，同时考虑了物质跃迁需克服一定能垒以及跃迁前、后的晶界能状态这两种因素，并且也能够模拟给出物质有可能通过能量起伏向高能态跃迁的情况，因此物理意义更为明确。

对所有格点遍历处理至到给定的模拟时间步后终止模拟，并给出相应的平均晶粒尺寸随模拟时间延长的变化趋势以及微观结构。

2 结果与讨论

2.1 特定模拟条件下的晶粒长大过程

图3中所示为G = 1、J = 1以及kT = 1.0时，平均晶粒半径随模拟时间延长的变化趋势及其拟合曲线。可见，随着模拟时间的延长，平均晶粒半径呈抛物线形式逐渐变长。其中，晶粒长大速率在模拟初期相对较快，并随模拟时间的延长而逐渐趋于平缓，与实际晶粒长大过程相符。对平均晶粒尺寸随模拟时间延长的变化趋势进行拟合，所得拟合曲线与平均晶粒半径数值变化之间吻合程度良好，并可得到拟合曲线的模拟长大指数为0.493，接近理论长大指数(0.5)。

图3 在G=1，J=1以及kT=1.0的条件下平均晶粒半径随模拟时间延长的变化趋势及拟合曲线Fig.3 Change tendency and fitted curve of average grain radius with extending simulation time under conditions of G=1,J=1, and kT=1.0

图4所示为晶粒长大微观结构随模拟时间延长的演变过程的模拟结果。可见，随着模拟时间的延长，小晶粒逐渐变小，直至在晶界处消失，而大晶粒则通过吞噬周围小晶粒逐渐长大，并具有5～7条平直晶界，最终引起晶粒总数量减少，晶粒平均尺寸增加。另外，也可以发现，晶界在向其曲率中心移动的过程中逐渐趋于平直化，有效避免了采用元胞自动机方法模拟晶粒长大过程时容易形成锯齿状不平整晶界的问题[18−19]。同时，三晶界交角逐渐转变成120˚，表明晶界处界面张力随着模拟时间的延长而逐渐达到平衡状态。另外，图5所示为铌镁酸铅陶瓷在1 200 ℃退火4 h后的断口显微结构，通过对比可以发现，模拟得到的晶界分布拓扑特征与实验观察结果较为一致。

图4 在 G = 1，J = 1和 kT = 1模拟条件下历经不同模拟时间后得到的模拟晶粒显微组织: (a) tMCS=120; (b) tMCS=240;(c) tMCS=360; (d) tMCS=480Fig.4 Development of simulation microstructure with extending simulation time under conditions of G = 1, J = 1 and kT = 1.0:(a) tMCS= 120; (b) tMCS=240; (c) tMCS=360; (d) tMCS=480

图6所示为经历不同模拟时间后的相对晶粒尺寸的分布情况。由图6可看出，不同模拟时间下所得到的相对晶粒尺寸均呈对数正态分布；进一步拟合得到的相对晶粒尺寸分布曲线也基本相同，即具有时间不变性。另外，也没有出现特别粗大的晶粒，最大晶粒尺寸约为平均晶粒尺寸的2.5倍。因此，可以认为本研究采用的改进MC方法可模拟得到较为理想的正常晶粒长大过程。

图5 铌镁酸铅陶瓷在1 200 ℃退火4 h后的断口显微结构Fig.5 Fracture microstructure of lead magnesium niobate after being annealed at 1 200 ℃ for 4 h

2.2 不同模拟温度下的晶粒长大过程

图7所示为G = 1和J = 1时，平均晶粒半径在不同kT随模拟时间延长的变化趋势和拟合曲线。由图7可见，随着模拟温度的升高，晶粒长大速度呈加快趋势，无需重新构建温度—时间对应关系就能够较好地修正原有MC模拟方法中出现晶粒长大速率随温度升高而减慢这一与实际晶粒长大规律不符的模拟现象。另外，在不同模拟温度下，平均晶粒尺寸也均随模拟时间的延长而呈抛物线形式增长，拟合后可得到模拟晶粒长大指数分布在 0.478～0.493的范围内，均与理论长大指数(0.5)较为接近。

图6 G = 1，J = 1以及kT = 1.0时经历不同模拟时间后的相对晶粒尺寸分布情况Fig.6 Relative grain size distribution after simulation for different times under conditions of G = 1, J = 1 and kT = 1.0

图7 G = 1，J = 1，kT = 1.0、1.5和2.0时平均晶粒半径随模拟时间延长的变化趋势和拟合曲线Fig.7 Change tendencies and fitted curves of average grain radius with extending simulation times under conditions of G = 1, J = 1, kT = 1.0 (a), 1.5 (b) and 2.0 (c)

图8 G = 1，kT =1.0、1.5、2.0时不同模拟时间下的微观结构Fig.8 Simulation of microstructures at different extending simulation times under condition of G = 1: (a) kT =1.0, tMCS=240;(b) kT =1.0, tMCS=480; (c) kT =1.5, tMCS=240; (d) kT =1.5, tMCS=480; (e) kT =2.0, tMCS=240; (f) kT =2.0, tMCS=480

2.3 不同激活能垒状态下的晶粒长大过程

图9所示为在J=1，kT=1时平均晶粒半径随模拟时间延长的变化趋势和拟合曲线。由图9可见，在模拟时间相同时，平均晶粒半径随着激活能垒(G)的增加而逐渐减小。对平均晶粒半径数据进行拟合，可得到不同激活能垒状态(G)下的模拟晶粒长大指数在0.472～0.493的范围内。

图9 平均晶粒半径随模拟时间延长的变化趋势和拟合曲线Fig.9 Change tendency and fitted curve of average grain radius with extending simulation times

图10 kT = 1时不同模拟时间下的微观结构的模拟结果Fig.10 Simulation of microstructure at different extending simulation times under condition of kT=1: (a) G=1, tMCS=240; (b) G=1,tMCS=480; (c) G=2, tMCS=240; (d) G=2, tMCS=480

图10所示为J=1，kT = 1时，晶粒长大微观结构随模拟时间延长的演变过程。从10可看出，在模拟温度和时间条件相同时，随着激活能垒状态(G)的增加，物质跃迁至相邻其他晶粒并引起晶界迁移也变得更困难，并最终导致晶粒长大速率减慢，在模拟所得微观结构上表现为晶粒数量相应增加，晶粒尺寸逐渐减小。

3 结论

1) 在各种模拟条件下得到的模拟晶粒长大指数在 0.472～0.493的范围内，接近理论晶粒长大指数(0.5)。在模拟所得的微观结构中，小晶粒随模拟时间的延长逐渐变小，直至在晶界处消失，大晶粒通过吞噬周围小晶粒逐渐长大，晶界向其曲率中心移动并趋于平直化，三晶界交角逐渐转变成 120˚。通过改进MC方法可以模拟得到与晶粒长大动力学相关理论和实验现象符合较好的晶粒形貌演变过程和晶界分布拓扑特征。

2) 在跃迁概率计算中引入跃迁能垒因素，可将物质在晶界处的跃迁概率与模拟温度直接关联，即跃迁概率随模拟温度的升高而增加，进而能够模拟得到晶粒长大速率随模拟温度的升高而加快的模拟结果。另外，通过增加跃迁能垒，使物质跃迁至相邻其他晶粒并引起晶界迁移的概率降低，因此可以给出晶粒长大速率相应减慢的模拟结果。改变模拟温度和能垒等模拟条件后，相应得到的晶粒长大速率变化趋势与实际晶粒长大规律较为一致。

REFERENCES

[1] MORON C, MORA M, GARCIA A. Computer simulation of grain growth kinetic[J]. Journal of Magnetism and Magnetic Materials, 2000, 215/216: 153−155.

紫趾综合征的确诊需作受累器官的组织学样本检查，可在直径为100μm～200μm的小动脉血管的管腔内见到胆固醇晶体，镜下可见胆固醇晶体溶解并出现针状裂缝[4]。

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

[3] 魏承炀, 高英俊, 张丽娜. 晶粒长大的 Monte Carlo模拟方法—递归统计法测定晶粒度[J]. 中国有色金属学报, 2008,18(1): 132−137.WEI Cheng-yang, GAO Ying-jun, ZHANG Li-na. Monte Carlo simulation of grain growth—Recursive statistics method of grain size[J]. The Chinese Journal of Nonferrous Metals, 2008, 18(1):132−137.

[4] HUANG C M, JOANNE C L, PATNAIK B S V,JAYAGANTHAN R. Monte Carlo simulation of grain growth in polycrystalline materials[J]. Applied Surface Science, 2006, 252:3397−4002.

[5] ZHOU Ming, LI Shi-chen, ZHENG Zi-qiao, YANG Pei-yong.Kinetics Monte Carlo simulation of microalloying effect in Al-Ag alloys[J]. Transactions of Nonferrous Metals Society of China, 2007, 17(3): 461−467.

[6] 张继祥, 关小军, 孙 胜. 一种改进的晶粒长大 Monte Carlo模拟方法[J]. 金属学报, 2004, 40(5): 457−461.ZHANG Ji-xiang, GUAN Xiao-jun, SUN Sheng. A modified Monte Carlo method in grain growth simulation[J]. Acta Metallurgica Sinica, 2004, 40(5): 457−461.

[7] YU Q, ESCHE S K. A Monte Carlo algorithm for single phase normal grain growth with improved accuracy and efficiency[J].Computational Materials Science, 2003, 27: 259−270.

[8] 叶日晴, 赵建华, 何陵辉. 退火过程中晶粒长大的二维计算机模拟[J]. 无机材料学报, 2001, 16(1): 122−128.YE Ri-qing, ZHAO Jian-hua, HE Ling-hui. Computer simulation of grain growth in two dimensions in annealing process[J].Journal of Inorganic Materials, 2001, 16(1): 122−128.

[9] 陈健美, 周卓夫, 唐建国, 张新明. 一种考虑晶界能各向异性模拟晶粒长大的Monte Carlo方法[J]. 材料导报, 2003, 17(8):77−79.CHEN Jian-mei, ZHOU Zhuo-fu, TANG Jian-guo, ZHANG Xin-ming. A Monte Carlo simulation of grain growth with grain boundary energy anisotropy considered[J]. Materials Review,2003, 17(8): 77−79.

[10] OGIBAYSHI S. Monte Carlo simulation of grain growth taking into account the influence temperature[J]. ISIJ International,2008, 48(3): 368−376.

[11] RADHAKRISHNAN B, ZACHARIA T. Monte Carlo simulation of grain boundary pinning in the weld heat-affected zone[J].Metallurgical and Materials Transactions A, 1995, 26(8):2123−2130.

[12] RADHAKRISHNAN B, ZACHARIA T. Simulation of curvature-driven grain growth by using a modified Monte Carlo algorithm[J]. Metallurgical and Materials Transactions A, 1995,26(1): 168−180.

[13] GAO J H, THOMPSON R G. Real time-temperature models for Monte Carlo simulations of normal grain growth[J]. Acta Materialia, 1996, 44: 4565−4570.

[14] GEIGER J, ROOSZ A, BARKOCZY P. Simulation of grain coarsening in two dimensions by cellular automaton[J]. Acta Materialia, 2001, 49: 623−629.

[15] BINDER K. Methods in statistical physics[M]. Heidelberg:Springer, 1986.

[16] 刘祖耀, 李世晨, 郑子樵, 陈大钦, 李 剑. 正常晶粒长大的计算机模拟(Ⅰ)—晶粒长大动力学跃迁概率的改进[J]. 中国有色金属学报, 2003, 13(6): 1357−1360.

LIU Zu-yao, LI Shi-chen, ZHENG Zi-qiao, CHEN Da-qin, LI Jian. Computer simulation of grain growth (Ⅰ)—Modified transition probability[J]. The Chinese Journal of Nonferrous Metals, 2003, 13(6): 1357−1360.

[17] 王海东, 张 海, 李海亮, 汤育才. 焙烧过程晶粒长大的

Monte Carlo模拟[J]. 中国有色金属学报, 2007, 17(6):990−996.

WANG Hai-dong, ZHANG Hai, LI Hai-liang, TANG Yu-cai.Monte Carlo simulation of grain growth in calcination process[J].The Chinese Journal of Nonferrous Metals, 2007, 17(6):990−996.

[18] 麻晓飞, 关小军, 刘云腾, 申孝民, 王丽君, 宋述同, 曾庆凯.基于改进转变规则的晶粒长大 CA模型[J]. 中国有色金属学报, 2008, 18 (1): 138−144.MA Xiao-fei, GUAN Xiao-jun, LIU Yun-teng, SHEN Xiao-min,WANG Li-jun, SONG Shu-tong, ZENG Qing-kai. Cellular automaton model for grain growth based on modified transition rule[J]. The Chinese Journal of Nonferrous Metals, 2008, 18(1):138−144.

[19] HE Y Z, DING H L, LIU L F, SHIN K. Computer simulation of 2D grain growth using a cellular automata model based on the lowest energy principle[J]. Mater Sci Eng A, 2006, 429:236−246.

(编辑 杨 华)

Modified Monte Carlo simulation of normal grain growth considering energy barrier for mass transition

LIU Jian-yuan, ZHANG Xiao-yong, ZHOU Ke-chao
(State Key laboratory of Powder Metallurgy, Central South University, Changsha 410083, China)

TG 111; TP 391.7

A

1004-0609(2010)05-0961-08

国家重点基础研究发展计划资助项目(2005CB623703)；中南大学博士后科学基金资助项目(2009年)

2009-07-15；

2009-12-17

张晓泳，博士；电话：0731-88830464；E-mail: zhangxiaoyong@mail.csu.edu.cn