晶粒生长的蒙特卡罗模拟

计算材料学课程设计指导老师:江建军 教授电子科学与技术系2004年6月晶粒生长的 Monte Carlo 模拟梁树雄，韩轲，张士亮，柏帆，胡斯杨，陈嵩，蒋雄军，严军刚，程晨（华中科技大学电子科学与技术系,湖北 武汉 430074）摘要：在实际生产中，晶粒的长大现象非常普遍的存在于金属、合金的凝固，陶瓷的烧结等 过程中，直接影响着所得材料的性能。所以研究晶粒生长是材料工程中的重要课题，但直接 研究速度慢，容易受到条件限制，而计算机模拟具有速度快，数据准确，成本低和直观简便 等一系列优点，经多年发展，已成为除实验和理论外的第三种关于晶体生长的研究手段，有 着不可比拟的优越性。 对于晶体生长来说， 所用的模型及方法大致分为蒙特卡罗方法 （Monte Carlo）简称MC方法或改进的MC方法。关键词：蒙特卡罗方法；晶体生长；计算机模拟Monte Carlo Method simulation of the grain growthAbstract : Grain growth exists very frequently in the condense and amalgamation of metal or alleys, processes of sintering of pottery. As its position, the capability of the materials attribute directly to the phenomena. With the great topic of materials engineering, however, due to restriction conditions, no more progress has achieved in the research. Nowadays, with the rapid progress of computational simulation and the improvement of the exactness. Besides of the traditional experiment and ordinary theoretic analogy, computational simulation has been regard as the other method to the problem. All of the models and methods can be divided to Monte Carlo Method （short for MC） and Improved MC Method.KEY WORDS ：Monte Carlo Method ；grain growth ；computational simulation一、引言晶粒凝聚的理论过程：高温的液态晶体在降温后即会凝聚成固体晶体。（1）成核现象 液态中随机在低温区凝聚的经理，发生在一个很小的过冷度范围内。（ 2）晶体生长液态晶体产生核化中心后，核化中心周围的液体开始逐渐凝固和生长。围绕着核化中心 即种子生长，也称作晶粒长大。晶粒长大可以分正常晶粒长大和不正常晶粒长大。正常晶粒 长大的基本特征是长大过程中晶粒尺寸保持基本均匀，相对晶粒尺寸分布具有不变性，通常 称之为自相似性。不正常晶粒长大是指晶粒长大由于某种原因，如第二相粒子的钉轧作用而 使长大受到限制，在某些条件如温差升到某临界值时，少数晶粒吞并其他晶粒而长大。（ 3）生长过程初期，在液态环境下围绕晶核向四周无限制的等速生长知道碰到另一个生长的晶粒，二 者中间形成晶界，此时停止生长。生长出的晶粒多不规则，晶界也没有特定的形态。后期， 晶粒以吞食的形式逐渐长大，以牺牲小晶粒为代价。随着时间的推移，晶粒的晶界不断变化，当时间足够长时，晶界多呈直线，晶界的边点多呈三叉状在实际生产中，晶粒的长大现象非常普遍的存在与金属，合金的凝固，陶瓷的烧结过程 中，直接影响着所得材料的性能。因而，对材料的晶粒长大过程中的成核及晶粒生长进行设 计，并用人工干预的方法，如控制炉温，加热时间，烧制方式及过程对材料进行性能上的改 善，使其可以达到使用要求，这十分具有实际意义。对晶粒生长的模拟，近年来突出的实现手段是：广泛利用了计算机系统来进行大量的计 算和图形界面。模拟计算机模拟具有速度快，数据准确，成本低和直观简便等一系列优点， 经多年发展，已成为除实验和理论外的第三种关于晶体生长的研究手段，有着不可比拟的优 越性。对于教学来说，用计算机程序对传统的实验过程进行图形化的模拟，已成为最普遍， 效果最好的模拟，设计方式。计算机模拟晶粒生长的方法和思路越来越多，所用的模型及模拟方法大体可分3种：MC方法或改进的(MC)方法、使用连续扩散界面场模型以及采用的其它模型；如Languerre模型，将细胞状的晶粒结构看作 Languerre棋盘形布局来处理。但最好的方法仍是MC法。在过去的几十年里，主要应用计算机模拟技术定量地理解晶粒生长过程，尤其是运用MC技术模拟等温条件下二维和三维的晶粒长大过程，并且从中获取从解析方程中难以得到的晶粒生 长的拓扑学信息。在本文里，我们重点讨论MC及其改进算法Q-state potts算法在模拟正常晶 粒生长中的作用，而对于异常晶粒生长，其物理模型构造及模拟过程存在一定的复杂性，在 此我们先不以予讨论。相关资料可以查找文献资料 2。并通过计算机构造一个 Voronoi网格 模型来给出模拟的过程及结果统计。二、基本模型1基本理论对多晶材料，在晶核形成之后，其晶粒生长现象分为两种不同类型：正常晶粒生长和异 常晶粒生长：(1) 正常晶粒生长时，晶粒尺寸呈一致增长的趋势，主要特点是晶粒尺寸F (R/ R )和拓扑 分布函数P (Nc)不随时间改变。其中,R为晶粒半径,R为平均晶粒半径，Nc是晶粒边数。 在等温条件下，正常晶粒的晶粒生长可以推导出下抛物线型动力学方程：R = Ktn(1)m m或 R - R(t=o)= Bt (2)在(1)式中，t为时间，K为生长率常数，n晶粒生长指数，通过(1)式对时间求导以及变换后，在Hillert的理论推导3中，理论预测稳态晶粒生长系数n应该等于0.5,但是实验和 模拟结果表明，在不忽略晶体拓扑结构影响的情况下，应当总是n<0. 5的这是由于实际材料和理想材料中存在的差异造成的。当时间 t很长时,若晶体形貌出现R»R° （R0为初始时晶 粒半径）,式（2）可以和式（1）等效，且m= 1/n，b为一常数，也就是说，Rm - Rm（t=o）与时 间成正比关系。（2）对于异常晶粒生长，则是指在重结晶的显微结构中，一些晶粒的尺寸迅速增加，最大尺寸的晶粒以比算术平均速率大得多的速率增长。异常晶粒生长的动力学方程可以表示为：X = 1 - exp - g 其中,X为二次重取向的晶粒的面积分数,g （ t）为与时间有关的函数。通常g （ t） = a t2 ,即修正方程：X = 1 - bexp (- a t p)其中,b、p为常数，且目前大部分p的实测值为1. 8± 0. 3左右图1 V oronoi网格一般认为晶粒生长变化的直接原因是驱动力的改变。 晶粒生长的驱动力主要来自总的晶界能的减少。但由 于生长的复杂性，迁移率、表面能、以及片材中的曲率、 应力等因素也将引起附加驱动力，从而引起晶粒异常 生长。通过上边的晶粒的生长动力学理论，我们在下 边对正常晶粒生长进行建模和模拟。2 基本模型（1）二维正常晶粒生长初始状态模拟的构造：首 先将晶体平面拓扑结构用网络格点表示出来，在实际 构造中，把晶体的晶格微观结构划分为不连续的晶格（正方形，三角型或六边型）我们参照了二维结构细胞结构构造方法，将微观结构划分成六边形结构，并映射到一个离散的Voronoi网格中,如图1。初始图形的生成方法采用细胞几何扩散法：在平面上随机撒下一些晶向不同的格点， 称为核点，核点取是一个从1到Q的值，表明该点的晶粒取向。在这样的模型中，我们可以通 过几何计算确定每个核点在等速扩散生长之后，与其他核点的交界，称之为晶界，晶界定义 为具有不同取向的两个晶粒的公共边，并且每个晶粒都是由离散化的晶格组成。这样，一个晶体的成核，以及晶粒均匀生长，直到晶粒间相互接触而停止生长，并形成晶界的一系列过程就完成了。下一步就是通过这样的Voronoi网格初始图形模拟正常晶粒的后期生长过程，我们将采用MC算法来完成图2六边形格点(2) 经典的蒙特卡罗(MC )算法模拟过程：晶粒生长变化的显微结构的演化可由改变每一晶格的取 向数(自旋数翻转来实现)，其中自旋数就是我们在前边撒点时 的核点取向数Q。在我们构造好的网格当中，随机抽取一个 格点，如图2,计算它与周围所有近邻格点的相互作用一般认 为这种相互作用为MNE = -J » _1)ij其中，s是格点i的晶向数，其数值为Q个取向中的一个，S是格点i的近邻格点的晶粒晶向 数。S Sis j为delta函数,M,N为格点数。其中，计算出来的E就是这个格点所具有的晶界能_ ?1 当 Sj=s 时爲=?0当Sj芦时由式(5)可见，求和包括所有近邻格点。其实在计算晶界能时理论上应该计算所有格点对 选取的格点的相互作用，但是对于实际应用，我们近似看成远邻格点对该格点的作用相等且 远小于近邻格点的作用能，因此我们认为只有近邻格点才对试验格点的局部能量有贡献。运用经典的MC算法(6 )，我们模拟判断其每个格点的晶向变化步骤如下：(1) 随机选取一个试验格点，其晶格取向为al，计算其格点晶界能E。(3) 进行随机取向试验：在其余 Q - 1种晶格取值向中任选其一，赋予该晶格，即晶格的数 字改变为a2 (a2却1)，得到改变后的晶界能E2。(4) 计算取向改变后相应的能量改变 E。通过由式计算晶粒晶向改变的几率 W，判断取向改变是否可以实现。?exp(-? E/KbT) ?E >0W?八B '(7)?1?E O 若 E< 0;则认为取向改变是允许的； 若 E> 0,则取向概率为exp (- E/KBT)，此时，在0, 1 间产生一个均匀分布的随机数n 如果exp (- E/KBT) > n，则取向改变允许,否则不允许,即只有使体系能量减少的转变才 是允许的。模拟中，我们将判断N取向实验的次数作为为一个时间单元，称蒙特卡罗时间步长(MonteCarlo Step， MCS), 其中 N 为格点总数,。从 MCS 到真实时间存在一个转换，即具有一个 简单的活化能因子 e , 此因子对应于原子的跃迁频率。3MC 模拟实际算法 由于上述经典的蒙特卡罗算法的计算机模拟时间相对较长，在通过实践验证后发现普通PC 无法有效实现，因此实际上在晶粒生长过程中我们采用的算法是基于 MC 算法的 Q-state potts 算法及其改进算法：Q-state potts 算法的实现过程如下：(1) 首先随机选取一个格点，晶粒取向为 al计算所选格点的六个最近的相邻格点中相同取 向数字也为a的格点数目N1 ;(2) 将这个选定格点的取向数字，随机地变为1Q中的另一个数b,再计算与这个新取向数 字b相同的近邻格点数N2 ;(3) 如果N2 >=N1 ,则将上述选取的晶格确定为新的取向 b;否则,恢复原来的取向a;(4) 随机选取另一个格点 ,重复以上步骤。(5) 若Voronoi网格总格点个数为n,则每完成n次随机抽取定为一个蒙特卡罗步(Monte Carlo Step)Q-state potts 算法是建立在蒙特卡罗模拟基础上的改进算法。在简化能量的计算量的同时 保留了蒙特卡罗算法随机性模拟的本质 ,但是也还存在一些问题：(1) 在对随机格点的选取上，Q-state potts要浪费大量的计算时间来进行判断，其实只有抽样 点落在晶体的边界处时，才可能产生晶粒生长。而晶体的边界在实际晶体中所占面积很小 ,大 部分判断都落在都是晶体的内部。(2) 另外经过模拟验证，晶粒取向变成所有可能的取向Q 中的一种的概率是极低的 , 因此 对晶粒取向变化可能值的确定是应该可以减少到一定程度的。通过上边两缺点 对Q-state potts算法作改进，可以得到另一种改进的Q-state potts算法：(1) 首先在晶界上选取一个格点，晶粒取向为a1,计算所选格点的六个最近的相邻格点中相 同取向数字也为 a 的格点数目 N1 ;(2) 随机地选取此晶粒最近邻另一个晶格的晶向 b (b工，计算与这个新取向数字b相同的 最近邻格点数 N2 ;(3) 如果N2 >=N1 ,则将上述选取的晶格确定为新的取向 b;否则保留原来的取向a;(4) 随机选取晶界另一个格点 ,重复以上步骤。（5）若晶界总格点个数为n,则每完成n次随机抽取定为一个蒙特卡罗步（Monte Carlo Step） 改进的Q-state potts算法的特点主要依据了两个观点：（1） 晶粒取向数字发生变化的区域为晶粒的边界，如果在找到边界格点后再进行判断，可以 大大缩短计算时间，并且此类忽略也是符合实际情况的。（2） “近邻”优于“远邻”理论：近邻格点比远邻格点具有更低的能量，并可以假设后远邻 格点都具有相同的较低能量。依据近邻优于远邻理论，可以将取向数字判断变为：只取最近 邻格点的取向数字中的一个，而忽略掉远邻格点的影响。利用改进的Q-state potts算法编写的程序，能够节省大量的计算时间，在较短时间内完成需要 大量时间进行随机抽取并判断变化的晶粒生长的模拟，可以使基于 MC算法的晶粒生长演示 在中低档个人PC上得以顺利运行。三、计算机的处理与作图1 初始图形生成过程初始图形将用Voronoi网格表示215X 255个晶格的区域，每个六边形代表一个晶格，此时代表晶格晶向的数字均置零。仿真开始时，利用计算机提供的随机函数产生 2n个随机数（XY），即在网格区域内产生n个晶核。如图3模拟中我们撒下了 150个随机点。然后以这些晶核为核心，每个点都以匀速向四周扩散，扩散后的每个点的晶相都改变成晶核的晶向， 直到晶粒间相互接触而停止生长，每两个晶粒间的交界即为晶界（本模拟实际通过求几何垂 直平分线确定晶界）。此时晶体的初始生长完成（如图 3）。图3初始图形生成2对晶体的生长进行MC方法的成长模拟（如图46）图 4 ( 50MCS)图 5( 150MCS)图 6 (400MCS)3 结论分析和处理利用一组统计数据来进行模拟结果的数值统计:当前时间（:MCS ）晶向个数01493525014636010014137215013339520012641725011944130011246935010848640010251545098536500945595509257160088597平均面积65084625我们通过上述数据绘出晶粒平均面积S-t图线(图7)和晶粒平均半径R-t图线(图8),m mm m其中R近似取R =、. S。利用上述R - R (t=0) = Bt，知道R - R (t=0)与时间成正比关系。其中 我们取经验值n = 0.47，m=1/n代入验证，通过matlab绘图并拟合曲线得到图9曲线，可以 看出和理论还是比较接近的，Rm - Rm(t=0) t曲线呈线性关系变化，结论：本次晶体生长的 蒙特卡罗蒙模拟和真实情况相差不大，与实验数据较符合。砂血如60040)400M圧旧£店弹硏ia0 1002003004006006007001QQ xn 300400500 GOO TOO图8晶粒平均半径R-t图线m m,图 9 R - R(t=0)- t 图线以上为运用蒙特卡罗(Monte Carlo)方法模拟二维正常晶粒生长的方法。然而根据现有模拟算法对模拟晶粒生长的过程而言，仅仅是一个比较理想和约束性条件较多的模拟算法。晶粒生长是一个复杂的过程，要对其进行完全且准确的模拟，还有待进一步的研究。参考材料：1 钟晓征, 陈伟元, 王豪才 等. 多晶材料晶体生长的 Monte Carlo 计算机模拟方法 I 模拟正常晶粒生长. 功能材料 , 1999, 30 (3)：232 2352 钟晓征,陈伟元,王豪才 等.多晶材料晶体生长的Monte Carlo计算机模拟方法II模拟异常晶粒生长. 功能材料 , 1999, 30 (3)：2362383 Hillert M. Acta. Metall,1965,13:227-2384 傅廷亮.基于C语言的二维细胞结构蒙特卡罗法仿真微电子与计算机1998, 3: 53 555 温俊芹,刘新田. 蒙特卡罗方法及其在晶粒生长模拟中的研究进展. WeldingTechnology,2000, 31(2): 12 莫春立，丁春辉，何若宏用Monte Carlo 方法模拟晶粒长大.沈阳工业学院学报2001,20(1): 61-667 陈伟元, 钟晓征, 谭锐, 王豪才. 多晶材料晶粒生长的计算机模拟研究 . 原子与分子物理学报, 2000, 17(2): 297-302