文章编号：1004-0609(2008)01-0132-06
晶粒长大的Monte Carlo模拟方法
——递归统计法测定晶粒度
魏承炀，高英俊，张丽娜
(广西大学 物理科学与工程技术学院，南宁 530004)
摘 要：应用Monte Carlo(MC)法模拟在周期性边界条件下的晶粒长大行为。利用MC法模拟时，晶界处格点的迁移引起晶粒的长大，根据这一主要特征提出一种精确快速的测定晶粒度的新方法—递归统计法，然后采用递归统计方法测量晶粒度。结果表明，递归统计法测得的晶粒度比截点法的更精确，而且测量精确度不受模型的格点类型以及晶粒的尺寸、形状等的影响，测量速度比其他统计方法要快。
关键词：递归统计；晶粒度； Monte Carlo方法；晶粒长大
中图分类号：TG 111.7；TP 391.9 文献标识码：A
Monte Carlo simulation of grain growth
—Recursive statistics method of grain size
WEI Cheng-yang, GAO Ying-jun, ZHANG Li-na
(School of Physical Science and Engineering, Guangxi University, Nanning 530004, China)
Abstract: Monte Carlo method was used to simulate grain growth process under periodic boundary condition. A new improved measurement method of grain size named recursive statistics method was introduced according to the fact that grain growth induced by displacement of single lattice around the grain boundary, and then the recursive statistics method was used to measure grain size. The results show that the grain size measured by the recursive statistics method is more accurate than the one measured by intercept method; the lattice types of method, sizes and shapes of grains cannot affect the measurement accuracy of the recursive statistics method; and the measuring speed of statistics method is faster than those of other statistics methods.
Key words: recursive statistics method; grain size; Monte Carlo method; grain growth
在纯金属、合金、陶瓷等多晶材料中，晶粒长大是最普遍的现象，对材料的性能有很重要的影响。目前应用最多且较成熟的模拟晶粒长大的方法是Monte Carlo(MC)法[1?4]，这是因为MC模拟可以揭示材料组织演变时拓扑学和动力学特征，可以进一步将其与实际的工艺过程结合，模拟在一定的工艺条件下的材料微观组织的演变情况。
在MC模拟晶粒长大、分析模拟结果的研究中，统计晶粒尺寸是一个必不可少的环节。目前，对晶粒数目的统计通常采用截点法[5]、快速标记算法[6]和图像跟踪技术。张继祥，关小军等[7]用近似定量金相方法即截点法计算出晶粒尺寸；莫春立等[6]利用光标在模型中按行列逐点标记和修正的方法即快速标记算法计算出晶粒尺寸。经过分析，截点法近似测量[5]模拟实验的晶粒度，没有充分发挥出计算机模拟实验的优势，而快速标记算法又因逐点修正引起了计算量的大幅度增大。因此，本文作者根据模拟时模型的格点迁移形成晶粒以及晶界这一主要特征[8?13]，提出了一种高效快速的精确统计算法，即递归统计法。这种算法利用统计模型中格点数，因而能精确统计出模拟过程的晶粒度，弥补了截点法在模拟实验中测量晶粒度的不足；并且采用递归游走的方式测量晶粒度，不需要修正计算，大大减小了计算量。
1 Monte Carlo模拟晶粒长大算法的描述
美国Exxon小组首先将Monte Carlo方法应用到晶粒长大过程模拟[1?4]，其主要内容如下。
将所选定的二维平面区域离散为若干微小的正多边形单元(如三角形、四边形或六边形)，对于每个单元，从Q(Q＞1)个整数中随机地选取一个作为其微观取向，相邻的相同取向的小单元构成一个晶粒，相邻的不同取向的单元之间形成晶界(如图1所示)。晶界能由选定的单元与其若干最相邻的单元的微观取向组合来定义，可描述为
晶界能的减小引起晶界的迁移。其算法的实现过程为
1) 确定单元的结构，对所计算的二维区域进行单元划分。
2) 将该区域映射到一个二维矩阵中，并对该矩阵进行初始化(给该矩阵的每一个元素随机地赋予一个1~Q的整数，该整数代表其对应单元的取向)。
3) 随机地选取一个单元，计算该单元的自由能E0。
4) 从其他Q~1个可能的取向中随机地选取一个赋给该单元并且计算其自由能E1。
5) 比较单元取向改变前后的能量差?E(?E=E1?E0)。新取向是否被接受，由下式判断：
N个这样的再定向尝试就构成了一个Monte Carlo步(MCS)。
6) 晶界处一系列单元取向的转变就构成了晶界的迁移，其迁移速度可由下式表示：
式(3)与经典晶粒长大速率是一致的。
图 1 二维矩阵单元的晶粒结构
Fig.1 Grain structure of 2-D matrix (The numerals express microcosmic direction of unit grid)
2 晶粒尺寸的精确快速统计
在模拟过程中，需要随时对晶粒的大小进行统计。但由于模拟过程中，晶粒是随机长大的，其尺寸、形状和分布无规律，有的呈现多边形，有的呈仿锤形，并且晶粒之间互相嵌入给晶粒统计带来了困难。为了方便统计，在用本文介绍的递归统计法统计晶粒尺寸之前需先对模型进行分级。
2.1 模型的分级
从MC模拟方法的基本思想和算法可以知道，晶粒长大时晶界处的格点迁移是该模拟方法的主要特 征[6?13]；晶粒长大过程都是由模型点阵内一系列单元取向的转变构成晶界的迁移而实现的；模型中任何形状的单位晶粒都是由一系列的点阵格点组成的。测定晶粒度，实质是统计各单位晶粒所占的格点数。递归统计法是基于模型格点(例如potts格点)的特征以及晶粒模拟长大时只考虑最近邻的格点的影响等特征提出的，因此，为了方便统计格点数，首先把单位晶粒占据的所有格点，按照被统计的先后顺序分为一级格点、二级格点、三级格点、直至n级格点。
模型的分级是为了能在统计时使统计步骤形成递归关系。本文以二维蜂窝状格点为例，介绍在晶粒统计计算时，对晶粒占据的所有格点进行分级的方法(分级方法如图2所示)和优点：
1) 统计每个单位晶粒时，第一个统计到的格点即为源格点，源格点周围的所有最近邻的格点为一级格点，一级格点最近邻除去已分级格点的所有格点称为二级格点，同理可对整个单位晶粒占据的所有格点进行分级；
2) 源格点最近邻最多有6个比自身低一级的格点，而一级及更低的格点最近邻最多只有3个比自身低一级的格点；
3) n级格点(也称末级格点)最近邻没有比自身低一级的格点；
4) 对模型分级后，统计计算晶粒度时可以改变传统的行列统计模式，而采用按格点的等级从高到低游走模式统计；
5) 由于除中心源格点外，其他格点的性质完全一样，统计步骤和条件也必然完全一样，因此该统计过程可以利用递归原理进行统计；
6) 分级时，容易对模型边界的格点引入周期性边界条件，其周期边界条件的引入方法与MC模拟过程相似。
总之，为了能够使模型的格点在统计时形成递归关系，必须对不同的格点类型进行适当的分级，分级后的格点形成图2所示的关系，按照其级别高低形成从小到大的包围圈。本文提出的递归统计方法正是应用各个包围圈的性质相同这一特点使统计过程形成递归关系。对于不同的格点类型，只需对格点按照前面介绍的方法进行分级后，就可以应用递归统计法进行晶粒统计；对于三维点阵模型，也可以用同样的方法进行分级，分级后仍然可以用同样的统计方法统计晶粒度，只是各个包围圈转变成立体包围。
图2 模型的分级
Fig.2 Classification of model
2.2 递归统计方法
对晶粒统计分析时，递归统计法是利用光标在模拟结果的数值阵点中游走时形成的递归关系，然后用计算机实现这种递归关系而统计出阵点数目的方法。
晶粒递归统计的具体步骤如下。
a. 首先把晶粒中的任意一个格点定为中心源格点，晶粒占据的所有格点就以该格点为基准进行分级。中心源格点选定后，晶粒占据的所有格点的性质也就确定下来。用传统的行列游走统计一个晶粒尺寸未知的晶粒时，由于晶粒的形状不规则，所以游走的退出条件不确定，就产生了统计过程中繁琐的修正，导致了计算量增大等缺点。然而递归统计法，首先对晶粒占据的所有格点进行分级，然后根据格点等级规定统计游走的每条路径都从高级的格点一直游走到最低级的格点才退出，也就是说光标只有到达晶粒的边界才会结束统计。
b. 以蜂窝状模型为例对统计过程进行说明。首先按照前文介绍的分级方法对模型格点进行分级。假设二级格点即是末级格点，统计过程的游走路径如下：
1) 从中心源格点开始游走到第一个一级格点；
2) 判断该格点是否是末级格点，如果是末级格点，则执行3)步，否则执行4)步；
3) 游走到下一个一级格点，并执行2)步；
4) 游走到第一个二级格点；
5) 判断该格点性质，如果是末级格点，则执行6)步，否则继续向低级格点游走；
6) 如果仍有二级格点没有统计，则游走到下一个二级格点，并执行5)步，否则执行3)步。
c. 统计过程形成递归关系。本文把统计过程2)、3)步的循环过程称为一级格点循环过程；第5)、6)步的循环过程称为二级格点循环过程。依此类推，有三级格点循环过程、四级格点循环过程……n级格点循环过程。
图3所示为步骤b。可以看出，二级格点循环过程是包含在一级格点循环过程中。在统计晶粒时，如果末级格点不是二级格点，那么这样包含着子循环的循环过程就不止一个，这些相互包含的循环过程形成了递归循环关系。这种递归循环关系将一直保持到末级格点才会消除，循环才会随之停止，这样才算是完成了一次递归过程。于是，统计过程就从中心源格点一直游走到晶粒的边界，并且对统计过的格点进行标记。结束一次递归过程后，重新从未标记过的格点中任选一个格点为源格点继续进行递归统计，直到全部格点都被标记为止，晶粒度的统计也才完成。
图3 递归统计流程图
Fig.3 Flowchart of recursion statistics
3 模拟结果分析
采用600×600 Potts (菠茨)格点。每一个格点都给定一个取向Q (Q的取值范围为1~64)，并引入周期性边界条件。图4所示为晶粒组织随时间的延续不断长大演变过程的模拟结果。可以看出，随着模拟时间的增加，平均晶粒度明显增大，晶粒长大是大晶粒吞噬小晶粒的结果；晶界交点处大都是三晶界相交，交角基本上等于120?。
为了得到晶粒长大指数n的精确数值，可以由 lg Rm—lg t直线的斜率求出[14]。图5(a)所示为本文模拟结果得到的lg Rm—lg t曲线，基本上接近于一条直线，其中，R的计算采用本文介绍的递归统计法，R的单位为μm，t的单位为MCS。根据最小二乘法将模拟数据耦合为直线，从直线的斜率得出平均晶粒长大指数n=0.460 3。
宋晓艳等[15]对MC模拟方法作了改进，改进后平均晶粒长大指数(≥0.485)比Anderson等[1]的模拟结果(0.414±0.03)更加接近于理论值(0.5)[16]。本文采用经宋晓艳等改进过的模拟方法，晶粒度的测量方法使 用本文介绍的递归统计法，测得的晶粒长大指数
图4 晶粒长大组织演变的模拟
Fig.4 Simulation images of microstructure evolution of grain growth: (a) 500 MCS; (b) 2 000 MCS; (c) 5 000 MCS; (d) 10 000 MCS
图 5 lg Rm—lg t的模拟曲线及其拟合直线
Fig.5 Simulation curves of lg Rm—lg t and its fitting line: (a) Simulation curve of lg Rm—lg t; (b) Simulation curve using recursion and intercept method
(0.460 3)与文献[15]中用传统的测定方法得到的模拟的结果 (≥0.485)有偏差。
图5(b)所示为分别用递归统计法与截点法测定晶粒度结果的比较。从图中可以看出，对于同一个模拟结果进行晶粒度的测量，用截点法测的数值要比递归统计法测得的结果大。因为模拟时间比较短时，晶粒的数目比较多，而且在所有晶粒中不规则的晶粒占了多数，这时采用截点法测量的晶粒数目就会比模拟的实际结果要偏多，所以得到的晶粒半径比模拟实际结果偏小；到了模拟后期，随着晶粒长大，晶粒数目越来越少，晶粒形状也变得较为规则，这时用截点法测得的结果，才较为可靠。从总的测量结果来看，用截点法得出的晶粒长大指数比实际模拟值要偏大，测量到n值为0.488 6，其测量到的n值要比模拟结果的实际数值要大得多，其偏差程度受晶粒形状的影响比较明显。而本文所采用的新计算方法是利用统计模拟结果中每个晶粒在点阵中所占据单元格的总个数而得到晶粒度，所以测量到的n值更接近于实际模拟值。
4 结论
1) 结合MC模拟晶粒长大时边界迁移的特征以及所应用的模型特点，提出了一种测量晶粒度的新方法——递归统计法。
2) 该方法比截点法能更精确地测量晶粒度，而且测量速度比截点法要快；不受模型点阵类型的限制，四方、六方和三角等点阵均适用。
2) 该方法对形状不规则的晶粒仍能得出准确的统计结果，并适合于二维和三维点阵模型；在有无周期性边界条件下均可精确计算晶粒个数、尺寸等参数。
REFERENCES
[1] ANDERSON M P, SROLOVITZ D J, GREST G S, SAHNI P S. Computer Simulation of Grain Growth—Ⅰ. Kinetics[J]. Acta Metal, 1984, 32(5): 783?791.
[2] SROLOVITZ D J, ANDERSON M P, SAHNI P S, GREST G S. Computer simulation of grain growth—Ⅱ. Grain size distribution, topology, and local dynamics[J]. Acta Metal, 1984, 32(5): 793?802.
[3] SROLOVITZ D J, ANDERSON M P, GREST G S, SAHNI P S. Computer simulation of grain growth—Ⅲ. Influence of a particle dispersion[J]. Acta Metal, 1984, 32(9): 1429?1438.
[4] GREST G S, SROLOVITZ D J, ANDERSON M P. Computer simulation of grain growth—Ⅳ. Anisotropic grain boundary energies[J]. Acta Metal, 1985, 33(3): 509?520.
[5] GB/T 6394—2002. 金属平均晶粒度测定方法[S].
GB/T 6394—2002. Metal-methods for estimating the average grain size[S].
[6] 张继祥, 关小军, 孙 胜. 一种改进的晶粒长大Monte Carlo模拟方法[J]. 金属学报, 2003, 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, 2003, 40(5): 457?461.
[7] 莫春立, 王维龙. 周期性边界条件下MC模拟晶粒数的统计[J]. 沈阳工业学院学报, 2003, 22(4): 36?39.
MO Chun-li, WANG Wei-long. Statistics of grain number in MC simulation under periodical boundary condition[J]. Journal of Shengyang Institute of Technology, 2003, 22(4): 36?39.
[8] 陈健美, 周卓夫, 唐建国, 张新明. 一种考虑晶界能各向异性模拟晶粒长大的Monte Carlo方法[J]. 材料导报, 2003, 17(8): 77?79, 60.
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, 60.
[9] 刘祖耀, 郑子樵, 陈大钦, 李世晨. 正常晶粒长大的计算机模拟机模拟(Ⅰ)—晶粒长大动力学跃迁概率的改进[J]. 中国有色金属学报, 2003, 13(6): 1357?1360.
LIU Zu-yao, ZHENG Zi-qiao, CHEN Da-qin, LI Shi-chen. Computer simulation of grain growth (Ⅰ)—Modified transition probability[J]. The Chinese Journal of Nonferrous Metals, 2003, 13(6): 1357?1360.
[10] 刘祖耀, 郑子樵, 陈大钦, 李世晨. 正常晶粒长大的计算机模拟机模拟(Ⅱ)—第二相粒子形状及取向的影响[J]. 中国有色金属报, 2004, 14(1): 122?126.
LIU Zu-yao, ZHENG Zi-qiao, CHEN Da-qin, LI Shi-chen. Computer simulation of grain growth (Ⅱ)—Influence of shape and orientation of second-phase particles on grain growth[J]. The Chinese Journal of Nonferrous Metals, 2004, 14(1): 122?126.
[11] 王 超, 刘国权. 晶粒长大过程中1种自相似态的仿真实验验证[J]. 稀有金属材料与工程, 2004, 33(2): 128?131.
WANG Chao, LIU Guo-quan. Simulation verification of self-similar state during process of normal grain growth[J]. Rare Metal Materials and Engineering, 2004, 33(2): 128?131.
[12] SOUCAIL M, M ESSINA R, COSNUAU A. Monte Carlo simulation of zener pinning in two dimensions[J]. Mater Sci Eng A, 1999, A271(1/2): 1?7.
[13] 张继祥, 关小军, 孙 胜, 申孝民, 董安平, 刘运腾. 晶粒长大过程微观组织演变Monte Carlo 方法模拟[J]. 山东大学学报(工学版), 2005, 35(4): 1?5.
ZHANG Ji-xiang, GUAN Xiao-jun, SUN Sheng, SHEN Xiao-min, DONG An-ping, LIU Yun-teng. Monte Carlo simulation of microstructure evolution during grain growth[J]. Journal of Shandong University (Engineering Science), 2005, 35(4): 1?5.
[14] MAO W M, ZHAO X B. Metal recrystallization and grain growth[M]. Beijing: Metallurgical Industry Press, 1994: 226?321.
[15] SONG X Y, LIU G Q. Modified Monte Carlo method for grain growth simulation[J]. Natural Science, 1998, 8(1): 92?97.
[16] CAHN J W, VLECK E V. Quadri junctions do not stop two-dimensional grain growth[J]. Scrip Mater, 1996, 34(6): 909?912.
基金项目：国家自然科学基金资助项目(50061001，50661001)；广西科学基金资助项目(桂科基0342004-1，桂科基0639004)
收稿日期：2007-06-28；修订日期：2007-11-20
通讯作者：高英俊，教授；电话：0771-3232666；E-mail: gaoyj@gxu.edu.cn
(编辑 陈爱华)
