材料设计与计算(一)

材料计算与设计绪论绪论1.1概述1.2材料设计的范围层次1.3材料设计的主要途径1.4材料设计面临任务1.1概述材料科学发展面临的问题1）由于研究对象的复杂性，现有理论手段很难处理一些极为复杂的问题，求解一个比较复杂的分子的薛定谔方程都很难实现；2）新的实验手段、仪器、设备虽然不断涌现，在一定范围内为实验研究提供了新方法，但大都极为昂贵，只为个别或少数拥有，研究的问题也极为有限。1.1概述材料设计计算的意义

1）计算机可以模拟进行现实中不能或很难实现的实验,如材料在极端压力、温度条件下的相变；

2）计算机可以模拟目前实验条件下无法进行的原子及以下尺度的研究等等；

3）计算机模拟可以验证已有理论和根据模拟结果修正或完善已有理论，也可以从模拟研究结果出发，指导、改善实验室实验。

因此,计算机模拟已成为除实验和理论外解决材料科学中实际问题的第3个重要组成部分,使材料的研究跳出了传统的“炒菜法”(trial-error)而发展为基于原理的方法。1.1概述1.1概述1.2材料设计的研究范畴按研究对象的空间尺度划分为三个层次:(1)工程设计层次:

尺度对应于宏观材料，涉及大块材料的加工和使用性能的设计研究。(2)连续模型尺度:

典型尺度在1μm量级，这时材料被看作连续介质，不考虑其中单个原子、分子的行为。(3)微观设计层次:

空间尺度在1nm量级,是原子、分子层次的设计"按空间和时间尺度分为四个层次:(1)宏观其空间尺度大致在0.1mm(目力能辨最小尺寸)至数万公里(人力跋涉之最远距离),时间尺度则大致在0.01秒(短跑时人所能分辨的速度最小差异)至100年(人的寿命差不多都在百年以内）(2)介观:

介观介于”宏观“与”微观“之间，其尺度主要在毫米量级。在材料学中其代表物是晶粒，也就是说需要注意微结构，如织构，成分偏析，晶界效应，孔中的吸附、逾渗、催化等问题。按空间和时间尺度分为四个层次:(3)微观其尺度主要在微米量级,也就是显微结构(世界)。以来借助于光学显微镜、电子显微镜、X衍射分析、电子探针等技术研究。在材料学中，这一尺度的代表物有晶须、雏晶、分相时产生的液滴等。(4)纳观:

其尺度范围在纳米至微米级，即10-6~10-9m，大致相当于几十个至几百个原子集合体的尺寸。1.3材料计算设计的主要途径材料数据库和知识库技术数据库一般包括材料的性能及一些主要的参数的数据：如材料成分、处理、试验条件以及材料的应用与评价。知识库主要是材料的成分、组织、工艺和性能间的关系及材料科学与工程的有关理论成果。它是实现人工智能的基本条件。实际上知识库就是材料计算设计中的一系列数理模型，用于定量计算或半定量描述的关系式。材料设计专家系统材料设计专家系统是指具有相当数量的与材料有关的各种背景知识。并能运用这些知识解决材料设计中有关问题的计算机程序系统。传统的专家系统主要有下列几个模块：优化模块、集成化模块、知识获取模块。现在逐步在发展智能专家网络系统是以模式识别和人工神经网络为基础的专家系统。材料系统设计材料科学将发展为材料系统科学,材料设计也必将是系统设计。计算机模拟技术它是利用计算机的计算推理功能，模拟实际情况来推测预报可能出现情况的一门技术。在80年代就进行了材料淬火过程的计算机模拟并建立了Metadex数据库1.4材料计算设计的任务(1)材料设计为国民经济和尖端技术服务，要结合国民经济建设和高技术项目开展材料设计工作。例如,厚壁压力容器材料、原子能应用材料、航空与航天用超高强度材料、高温合金!低温材料、电子信息材料、各种特殊功能材料等(2)进行多层次结合的材料计算设计。分层次研究的弱点是不同学科互相分割，难以取得系统的综合效果，特别是微观层次(电子、原子)的设计离预报、设计实际材料还有很大的距离。国内外很多有识人士提出了多层次综合的设计思想。(3)多学科的交叉、融合是必然的趋势。材料计算设计科学是材料、物理、化学、数学和计算机等多学科的交叉研究领域。鼓励材料科学和系统科学结合、整体化已成为当今科技发展的重要趋势，多层次和跨学科正是计算材料学的特点和本质。(4)数理模型的建立和实用化是关键。材料设计系统主要依赖于数理模型。各层次研究的关键是根据基础数据能否发展出符号实际的解析与数理模型，解决不同层次间计算方法的选择与整合。(5)材料计算设计科学的基础研究必须加强。我国的材料计算设计研究落后于国外,并且在观念!思维上也没有跳出国外现有的思路,有的还比较偏激。有关基础理论的研究工作有待于加强。第一章蒙特卡罗方法蒙特卡罗方法的基本思想蒙特卡罗方法的特点蒙特卡罗方法的主要应用范围蒙特卡罗方法的收敛性，误差作业第一节

蒙特卡罗方法概述第一章蒙特卡罗方法概述概念

蒙特卡罗方法：又称随机抽样方法或统计试验方法，是通过不断产生的随机序列来模拟过程。特点

1）蒙特卡罗方法是一种计算方法，但与一般数值计算方法有很大区别。它是以概率统计理论为基础的一种方法。

2）由于蒙特卡罗方法能够比较逼真地描述事物的特点及物理实验过程，解决一些数值方法难以解决的问题，因而该方法的应用领域日趋广泛。蒙特卡罗方法的基本思想

定义蒙特卡洛方法就是根据待求随机问题的变化规律，根据物理现象本身的统计规律，或者人为的构造一个合适的概率模型，依照大量的统计实验，使它的某些统计参量正好是待求解问题的解。

两个例子

例1.蒲丰氏问题

例2.射击问题（打靶游戏）基本思想计算机模拟试验过程例1.蒲丰氏问题为了求得圆周率π值，在十九世纪后期，有很多人作了这样的试验：将长为2l的一根针任意投到地面上，用针与一组相间距离为2a（

l＜a）的平行线相交的频率代替概率P，再利用准确的关系式：求出π值

其中Ｎ为投计次数，n为针与平行线相交次数。这就是古典概率论中著名的蒲丰氏问题。

一些人进行了实验，其结果列于下表：实验者年份投计次数π的实验值沃尔弗(Wolf)185050003.1596斯密思(Smith)185532043.1553福克斯(Fox)189411203.1419拉查里尼(Lazzarini)190134083.1415929例2.射击问题（打靶游戏）

设r表示射击运动员的弹着点到靶心的距离，ｇ(r)表示击中r处相应的得分数（环数），f(r)为该运动员的弹着点的分布密度函数，它反映运动员的射击水平。该运动员的射击成绩为

用概率语言来说，<g>是随机变量ｇ(r)的数学期望，即

现假设该运动员进行了Ｎ次射击，每次射击的弹着点依次为r1，r2，…，rN，则Ｎ次得分g(r1)，g(r2)，…，g(rN)的算术平均值代表了该运动员的成绩。换言之，为积分<g>的估计值，或近似值。在该例中，用Ｎ次试验所得成绩的算术平均值作为数学期望<g>的估计值（积分近似值）。

基本思想

由以上两个例子可以看出，当所求问题的解是某个事件的概率，或者是某个随机变量的数学期望，或者是与概率、数学期望有关的量时，通过某种试验的方法，得出该事件发生的频率，或者该随机变量若干个具体观察值的算术平均值，通过它得到问题的解。这就是蒙特卡罗方法的基本思想。当随机变量的取值仅为1或0时，它的数学期望就是某个事件的概率。或者说，某种事件的概率也是随机变量（仅取值为1或0）的数学期望。

因此，蒙特卡罗方法是用随机试验的方法计算积分，即将所要计算的积分看作服从某种分布密度函数f(r)的随机变量ｇ(r)的数学期望

通过某种试验，得到Ｎ个观察值r1，r2，…，rN（用概率语言来说，从分布密度函数f(r)中抽取Ｎ个子样r1，r2，…，rN，），将相应的Ｎ个随机变量的值g(r1)，g(r2)，…，g(rN)的算术平均值作为积分的估计值（近似值）。

存在问题：

1）为了得到具有一定精确度的近似解，所需试验的次数是很多的，通过人工方法作大量的试验相当困难，甚至是不可能的。2）本世纪四十年代以来，由于电子计算机的出现，使得人们可以通过电子计算机来模拟随机试验过程，把巨大数目的随机试验交由计算机完成，使得蒙特卡罗方法得以广泛地应用，在现代化的科学技术中发挥应有的作用。

计算机模拟试验过程

计算机模拟试验过程，就是将试验过程（如投针，射击）化为数学问题，在计算机上实现。以上述两个问题为例，分别加以说明。例1.蒲丰氏问题例2.射击问题（打靶游戏）由上面两个例题看出，蒙特卡罗方法常以一个“概率模型”为基础，按照它所描述的过程，使用由已知分布抽样的方法，得到部分试验结果的观察值，求得问题的近似解。例１．蒲丰氏问题

设针投到地面上的位置可以用一组参数（x,θ）来描述，x为针中心的坐标，θ为针与平行线的夹角，如图所示。任意投针，就是意味着x与θ都是任意取的，但x的范围限于［0，a］，夹角θ的范围限于［0，π］。在此情况下，针与平行线相交的数学条件是针在平行线间的位置

如何产生任意的（x,θ）？x在［0，a］上任意取值，表示x在［0，a］上是均匀分布的，其分布密度函数为：类似地，θ的分布密度函数为：因此，产生任意的（x,θ）的过程就变成了由f1(x)抽样x及由f2(θ)抽样θ的过程了。由此得到：其中ξ1，ξ2均为（0,1）上均匀分布的随机变量。每次投针试验，实际上变成在计算机上从两个均匀分布的随机变量中抽样得到（x,θ），然后定义描述针与平行线相交状况的随机变量s(x,θ)，为如果投针Ｎ次，则上式是针与平行线相交概率Ｐ的估计值。事实上，于是有例２．射击问题

设射击运动员的弹着点分布为用计算机作随机试验（射击）的方法为，选取一个随机数ξ，按右边所列方法判断得到成绩。这样，就进行了一次随机试验（射击），得到了一次成绩ｇ(r)，作Ｎ次试验后，得到该运动员射击成绩的近似值环数78910概率0.10.10.30.5蒙特卡罗方法的收敛性，误差

蒙特卡罗方法作为一种计算方法，其收敛性与误差是普遍关心的一个重要问题。收敛性误差减小方差的各种技巧效率收敛性

由前面介绍可知，蒙特卡罗方法是由随机变量X的简单子样X1，X2，…，XN的算术平均值：作为所求解的近似值。由大数定律可知，如X1，X2，…，XN独立同分布，且具有有限期望值（E(X)<∞），则即随机变量X的简单子样的算术平均值，当子样数Ｎ充分大时，以概率1收敛于它的期望值E(X)。误差

蒙特卡罗方法的近似值与真值的误差问题，概率论的中心极限定理给出了答案。该定理指出，如果随机变量序列X1，X2，…，XN独立同分布，且具有有限非零的方差σ2

，即

f(X)是X的分布密度函数。则

当N充分大时，有如下的近似式其中α称为置信度，1－α称为置信水平。这表明，不等式近似地以概率

1－α成立，且误差收敛速度的阶为。通常，蒙特卡罗方法的误差ε定义为上式中与置信度α是一一对应的，根据问题的要求确定出置信水平后，查标准正态分布表，就可以确定出。

下面给出几个常用的α与的数值：

关于蒙特卡罗方法的误差需说明两点：第一，蒙特卡罗方法的误差为概率误差，这与其它数值计算方法是有区别的。第二，误差中的均方差σ是未知的，必须使用其估计值来代替，在计算所求量的同时，可计算出。α0.50.050.003

0.67451.963减小方差的各种技巧

显然，当给定置信度α后，误差ε由σ和N决定。要减小ε，或者是增大N，或者是减小方差σ2。在σ固定的情况下，要把精度提高一个数量级，试验次数N需增加两个数量级。因此，单纯增大N不是一个有效的办法。另一方面，如能减小估计的均方差σ，比如降低一半，那误差就减小一半，这相当于N增大四倍的效果。因此降低方差的各种技巧，引起了人们的普遍注意。效率

一般来说，降低方差的技巧，往往会使观察一个子样的时间增加。在固定时间内，使观察的样本数减少。所以，一种方法的优劣，需要由方差和观察一个子样的费用（使用计算机的时间）两者来衡量。这就是蒙特卡罗方法中效率的概念。它定义为

，其中c是观察一个子样的平均费用。显然

越小，方法越有效。蒙特卡罗方法的特点优点能够比较逼真地描述具有随机性质的事物的特点及物理实验过程。受几何条件限制小。收敛速度与问题的维数无关。具有同时计算多个方案与多个未知量的能力。误差容易确定。程序结构简单，易于实现。缺点收敛速度慢。误差具有概率性。在粒子输运问题中，计算结果与系统大小有关。能够比较逼真地描述具有随机性质的事物的特点及物理实验过程从这个意义上讲，蒙特卡罗方法可以部分代替物理实验，甚至可以得到物理实验难以得到的结果。用蒙特卡罗方法解决实际问题，可以直接从实际问题本身出发，而不从方程或数学表达式出发。它有直观、形象的特点。受几何条件限制小

在计算s维空间中的任一区域Ds上的积分时无论区域Ds的形状多么特殊，只要能给出描述Ds的几何特征的条件，就可以从Ds中均匀产生N个点，得到积分的近似值。其中Ds为区域Ds的体积。这是数值方法难以作到的。另外，在具有随机性质的问题中，如考虑的系统形状很复杂，难以用一般数值方法求解，而使用蒙特卡罗方法，不会有原则上的困难。收敛速度与问题的维数无关

由误差定义可知，在给定置信水平情况下，蒙特卡罗方法的收敛速度为，与问题本身的维数无关。维数的变化，只引起抽样时间及估计量计算时间的变化，不影响误差。也就是说，使用蒙特卡罗方法时，抽取的子样总数N与维数s无关。维数的增加，除了增加相应的计算量外，不影响问题的误差。这一特点，决定了蒙特卡罗方法对多维问题的适应性。而一般数值方法，比如计算定积分时，计算时间随维数的幂次方而增加，而且，由于分点数与维数的幂次方成正比，需占用相当数量的计算机内存，这些都是一般数值方法计算高维积分时难以克服的问题。具有同时计算多个方案与多个未知量的能力

对于那些需要计算多个方案的问题，使用蒙特卡罗方法有时不需要像常规方法那样逐个计算，而可以同时计算所有的方案，其全部计算量几乎与计算一个方案的计算量相当。例如，对于屏蔽层为均匀介质的平板几何，要计算若干种厚度的穿透概率时，只需计算最厚的一种情况，其他厚度的穿透概率在计算最厚一种情况时稍加处理便可同时得到。另外，使用蒙特卡罗方法还可以同时得到若干个所求量。例如，在模拟粒子过程中，可以同时得到不同区域的通量、能谱、角分布等，而不像常规方法那样，需要逐一计算所求量。误差容易确定

对于一般计算方法，要给出计算结果与真值的误差并不是一件容易的事情，而蒙特卡罗方法则不然。根据蒙特卡罗方法的误差公式，可以在计算所求量的同时计算出误差。对干很复杂的蒙特卡罗方法计算问题，也是容易确定的。一般计算方法常存在着有效位数损失问题，而要解决这一问题有时相当困难，蒙特卡罗方法则不存在这一问题。程序结构简单，易于实现在计算机上进行蒙特卡罗方法计算时，程序结构简单，分块性强，易于实现。收敛速度慢

如前所述，蒙特卡罗方法的收敛速度为，一般不容易得到精确度较高的近似结果。对于维数少（三维以下）的问题，不如其他方法好。误差具有概率性由于蒙特卡罗方法的误差是在一定置信水平下估计的，所以它的误差具有概率性，而不是一般意义下的误差。在粒子输运问题中，计算结果与系统大小有关

经验表明，只有当系统的大小与粒子的平均自由程可以相比较时（一般在十个平均自由程左右），蒙特卡罗方法计算的结果较为满意。但对于大系统或小概率事件的计算问题，计算结果往往比真值偏低。而对于大系统，数值方法则是适用的。因此，在使用蒙特卡罗方法时，可以考虑把蒙特卡罗方法与解析（或数值）方法相结合，取长补短，既能解决解析（或数值）方法难以解决的问题，也可以解决单纯使用蒙特卡罗方法难以解决的问题。这样，可以发挥蒙特卡罗方法的特长，使其应用范围更加广泛。蒙特卡罗方法的主要应用范围

蒙特卡罗方法所特有的优点，使得它的应用范围越来越广。它的主要应用范围包括：粒子输运问题，统计物理，典型数学问题，真空技术，激光技术以及医学，生物，探矿等方面。随着科学技术的发展，其应用范围将更加广泛。蒙特卡罗方法在粒子输运问题中的应用范围主要包括：实验核物理，反应堆物理，高能物理等方面。蒙特卡罗方法在实验核物理中的应用范围主要包括：通量及反应率，中子探测效率，光子探测效率，光子能量沉积谱及响应函数，气体正比计数管反冲质子谱，多次散射与通量衰减修正等方面。作业

用蒲丰投针法在计算机上计算π值，取a=4、l=3。分别用理论计算和计算机模拟计算，求连续掷两颗骰子，点数之和大于6且第一次掷出的点数大于第二次掷出点数的概率。第二节随机数随机数的定义及产生方法伪随机数产生伪随机数的乘同余方法产生伪随机数的乘加同余方法产生伪随机数的其他方法伪随机数序列的均匀性和独立性作业第二节随机数

由具有已知分布的总体中抽取简单子样，在蒙特卡罗方法中占有非常重要的地位。总体和子样的关系，属于一般和个别的关系，或者说属于共性和个性的关系。由具有已知分布的总体中产生简单子样，就是由简单子样中若干个性近似地反映总体的共性。随机数是实现由已知分布抽样的基本量，在由已知分布的抽样过程中，将随机数作为已知量，用适当的数学方法可以由它产生具有任意已知分布的简单子样。随机数的定义及产生方法随机数的定义及性质随机数表物理方法随机数的定义及性质

在连续型随机变量的分布中，最简单而且最基本的分布是单位均匀分布。由该分布抽取的简单子样称，随机数序列，其中每一个体称为随机数。单位均匀分布也称为［0，1］上的均匀分布，其分布密度函数为：分布函数为：

由于随机数在蒙特卡罗方法中占有极其重要的位置，我们用专门的符号ξ表示。由随机数序列的定义可知，ξ1，ξ2，…是相互独立且具有相同单位均匀分布的随机数序列。也就是说，独立性、均匀性是随机数必备的两个特点。随机数具有非常重要的性质：对于任意自然数s，由s个随机数组成的s维空间上的点（ξn+1，ξn+2，…ξn+s）在s维空间的单位立方体Gs上均匀分布，即对任意的ai，如下等式成立：

其中P(·)表示事件·发生的概率。反之，如果随机变量序列ξ1,ξ2…对于任意自然数s，由s个元素所组成的s维空间上的点（ξn+1，…ξn+s）在Gs上均匀分布，则它们是随机数序列。由于随机数在蒙特卡罗方法中所处的特殊地位，它们虽然也属于由具有已知分布的总体中产生简单子样的问题，但就产生方法而言，却有着本质上的差别。随机数表

为了产生随机数，可以使用随机数表。随机数表是由0，1，…，9十个数字组成，每个数字以0.1的等概率出现，数字之间相互独立。这些数字序列叫作随机数字序列。如果要得到n位有效数字的随机数，只需将表中每n个相邻的随机数字合并在一起，且在最高位的前边加上小数点即可。例如，某随机数表的第一行数字为7634258910…，要想得到三位有效数字的随机数依次为0.763，0.425，0.891。因为随机数表需在计算机中占有很大内存，而且也难以满足蒙特卡罗方法对随机数需要量非常大的要求，因此，该方法不适于在计算机上使用。物理方法

用物理方法产生随机数的基本原理是：利用某些物理现象，在计算机上增加些特殊设备，可以在计算机上直接产生随机数。这些特殊设备称为随机数发生器。用来作为随机数发生器的物理源主要有两种：一种是根据放射性物质的放射性，另一种是利用计算机的固有噪声。一般情况下，任意一个随机数在计算机内总是用二进制的数表示的：其中εi（i=1，2，…，m）或者为0，或者为1。

因此，利用物理方法在计算机上产生随机数，就是要产生只取0或1的随机数字序列，数字之间相互独立，每个数字取0或1的概率均为0.5。用物理方法产生的随机数序列无法重复实现，不能进行程序复算，给验证结果带来很大困难。而且，需要增加随机数发生器和电路联系等附加设备，费用昂贵。因此，该方法也不适合在计算机上使用。伪随机数伪随机数伪随机数存在的两个问题伪随机数的周期和最大容量伪随机数在计算机上产生随机数最实用、最常见的方法是数学方法，即用如下递推公式：产生随机数序列。对于给定的初始值ξ1,ξ2…，ξk，确定ξn+k，ｎ=1，2，…。经常使用的是k=1的情况，其递推公式为：

对于给定的初始值ξ1，确定ξn+1，ｎ=１,２…伪随机数存在的两个问题用数学方法产生的随机数，存在两个问题：递推公式和初始值ξ1,ξ2…，ξk确定后，整个随机数序列便被唯一确定。不满足随机数相互独立的要求。由于随机数序列是由递推公式确定的，而在计算机上所能表示的[0，1]上的数又是有限的，因此，这种方法产生的随机数序列就不可能不出现无限重复。一旦出现这样的n＇，n″(n＇<n″)，使得下面等式成立：随机数序列便出现了周期性的循环现象。对于k=1的情况，只要有一个随机数重复，其后面的随机数全部重复，这与随机数的要求是不相符的。

由于这两个问题的存在，常称用数学方法产生的随机数为伪随机数。对于以上存在的两个问题，作如下具体分析。关于第一个问题，不能从本质上加以改变，但只要递推公式选得比较好，随机数间的相互独立性是可以近似满足的。至于第二个问题，则不是本质的。因为用蒙特卡罗方法解任何具体问题时，所使用的随机数的个数总是有限的，只要所用随机数的个数不超过伪随机数序列出现循环现象时的长度就可以了。用数学方法产生的伪随机数容易在计算机上得到，可以进行复算，而且不受计算机型号的限制。因此，这种方法虽然存在着一些问题，但仍然被广泛地在计算机上使用，是在计算机上产生伪随机数的主要方法。伪随机数的周期和最大容量

发生周期性循环现象的伪随机数的个数称为伪随机数的周期。对于前面介绍的情况，伪随机数的周期为n″－n＇。从伪随机数序列的初始值开始，到出现循环现象为止，所产生的伪随机数的个数称为伪随机数的最大容量。前面的例子中，伪随机数的最大容量为n″。产生伪随机数的乘同余方法

乘同余方法是由Lehmer在1951年提出来的，它的一般形式是：对于任一初始值x1，伪随机数序列由下面递推公式确定：

其中a为常数。乘同余方法的最大容量的上界对于任意正整数M，根据数论中的标准分解定理，总可以分解成如下形式：

其中P0=2，P1，…Pr表示不同的奇素数，α0表示非负整数，α1，…，αr表示正整数。a无论取什么值，乘同余方法的最大容量的上界为：

的最小公倍数。其中：关于a与x1的取值

如果a与x1满足如下条件：

对于

，

x1与M互素，则乘同余方法产生的伪随机数序列的最大容量达到最大可能值λ(M)。乘同余方法在计算机上的使用

为了便于在计算机上使用，通常取： Ｍ=2s

其中s为计算机中二进制数的最大可能有效位数

x1=奇数

a=52k+1

其中k为使52k+1在计算机上所能容纳的最大整数，即a为计算机上所能容纳的5的最大奇次幂。一般地，s=32时，a=513；s=48，a=515等。伪随机数序列的最大容量λ(M)=2s-2。

乘同余方法是使用的最多、最广的方法，在计算机上被广泛地使用。产生伪随机数的乘加同余方法产生伪随机数的乘加同余方法是由Rotenberg于1960年提出来的，由于这个方法有很多优点，已成为仅次于乘同余方法产生伪随机数的另一主要方法。

乘加同余方法的一般形式是，对任意初始值x1，伪随机数序列由下面递推公式确定：其中a和c为常数。乘加同余方法的最大容量关于乘加同余方法的最大容量问题，有如下结论：如果对于正整数M的所有素数因子P，下式均成立：当M为4的倍数时，还有下式成立：c与M互素，则乘加同余方法所产生的伪随机数序列的最大容量达到最大可能值M。

M，x1，a，c的取值为了便于在计算机上使用，通常取

M=2s

其中s为计算机中二进制数的最大可能有效位数。

a=2b+1 (b≥2) c=1

这样在计算中可以使用移位和指令加法，提高计算速度。

产生伪随机数的其他方法取中方法加同余方法伪随机数序列的均匀性和独立性

判断伪随机数序列是否满足均匀和相互独立的要求，要靠统计检验的方法实现。对于伪随机数的统计检验，一般包括两大类：均匀性检验和独立性检验。六十年代初，人们开始用定性的方法研究伪随机数序列的均匀性和独立性问题，简要叙述如下。伪随机数的均匀性

这里只考虑伪随机数序列ξ1,ξ2…，ξn全体作为子样时的均匀性问题。其中n为伪随机数序列的最大容量。对于任意的0≤x≤1，令Nn(x)表示伪随机数序列ξ1,ξ2…，ξn中适合不等式

ξi<xi=1，2，…，n

的个数，则标志伪随机数序列ξ1,ξ2…，ξn的均匀程度，称为均匀偏度。

将伪随机数序列ξ1,ξ2…，ξn从小至大重新排列并令，则由δ(n)的定义，容易证明很明显，对于固定的ｎ，δ(n)的值越小越好。它是描述伪随机数序列均匀程度的基本量。对于任意随机数序列，均有如下不等式成立：当时，所对应的伪随机数序列为最佳分布。

可以证明，伪随机数序列为最佳分布的充要条件是它取遍序列的所有值。对于计算机上使用的乘同余方法，按照前面介绍的方法选取a、x1时，所产生的伪随机数序列的均匀偏度对于乘加同余方法对于部分伪随机数的均匀性问题通常用统计检验方法检验。伪随机数的独立性

对于任意，令表示(ξ1,ξ2)，(ξ2,ξ3)，…，(ξn,ξn+1)中适合不等式

的个数，根据随机变量间相互独立的定义和频率近似概率的方法，令则ε(n)标志伪随机数序列ξ1,ξ2…，ξn的独立程度，简称为独立偏度。对于固定的n,ε(n)的值越接近于零，伪随机数序列的独立性越好。对于乘同余方法，对于乘加同余方法，因此，这两种方法的独立性都是很好的。同伪随机数的均匀性问题一样，伪随机数序列的独立性问题也是对它的全体讨论的。若只考虑伪随机数的一部分，在通常情况下给出ε(i)是相当因难的。因此，伪随机数序列的独立性问题的统计检验方法同样是非常重要的。作业

证明1—ξ是随机数。证明与同分布。第三节由已知分布的随机抽样随机抽样及其特点直接抽样方法挑选抽样方法复合抽样方法复合挑选抽样方法替换抽样方法随机抽样的一般方法随机抽样的其它方法作业第三节由已知分布的随机抽样

本章叙述由己知分布抽样的各主要方法，并给出在粒子输运问题中经常用到的具体实例。随机抽样及其特点

由巳知分布的随机抽样指的是由己知分布的总体中抽取简单子样。随机数序列是由单位均匀分布的总体中抽取的简单子样，属于一种特殊的由已知分布的随机抽样问题。本章所叙述的由任意已知分布中抽取简单子样，是在假设随机数为已知量的前提下，使用严格的数学方法产生的。为方便起见，用XF表示由己知分布F(x)中产生的简单子样的个体。对于连续型分布，常用分布密度函数f(x)表示总体的己知分布，用Xf表示由己知分布密度函数f(x)产生的简单子样的个体。另外，在抽样过程中用到的伪随机数均称随机数。直接抽样方法

对于任意给定的分布函数F(x)，直接抽样方法如下：其中，ξ1，ξ2，…，ξN为随机数序列。为方便起见，将上式简化为：若不加特殊说明，今后将总用这种类似的简化形式表示，ξ总表示随机数。证明

下面证明用前面介绍的方法所确定的随机变量序列X1，X2，…，XN具有相同分布F(x)。对于任意的n成立，因此随机变量序列X1，X2，…，XN具有相同分布F(x)。另外，由于随机数序列ξ1，ξ2，…，ξN是相互独立的，而直接抽样公式所确定的函数是波雷尔（Borel）可测的，因此，由它所确定的X1，X2，…，XN也是相互独立的（[P.R.Halmos,Measuretheory,N.Y.VonNosrtand,1950]§45定理2）。离散型分布的直接抽样方法

对于任意离散型分布：其中x1，x2，…为离散型分布函数的跳跃点，P1，P2，…为相应的概率，根据前述直接抽样法，有离散型分布的直接抽样方法如下：该结果表明，为了实现由任意离散型分布的随机抽样，直接抽样方法是非常理想的。例1.二项分布的抽样

二项分布为离散型分布，其概率函数为：其中，P为概率。对该分布的直接抽样方法如下：例2.泊松(Possion)分布的抽样

泊松(Possion)分布为离散型分布，其概率函数为：其中，λ>0。对该分布的直接抽样方法如下：例3.掷骰子点数的抽样

掷骰子点数X=n的概率为：选取随机数ξ，如则在等概率的情况下，可使用如下更简单的方法：其中［］表示取整数。例4.碰撞核种类的确定

中子或光子在介质中发生碰撞时，如介质是由多种元素组成，需要确定碰撞核的种类。假定介质中每种核的宏观总截面分别为Σ1，Σ2，…，Σn，则中子或光子与每种核碰撞的概率分别为：其中Σt＝Σ1＋Σ2＋…＋Σn。碰撞核种类的确定方法为：产生一个随机数ξ，如果则中子或光子与第I种核发生碰撞。例5.中子与核的反应类型的确定

假设中子与核的反应类型有如下几种:弹性散射，非弹性散射，裂变，吸收，相应的反应截面分别为Σel，Σin，Σf，Σa。则发生每一种反应类型的概率依次为：

其中反应总截面Σt＝Σel＋Σin＋Σf＋Σa。

反应类型的确定方法为：产生一个随机数ξ

连续型分布的直接抽样方法

对于连续型分布，如果分布函数F(x)的反函数

F－1(x)存在，则直接抽样方法是：例6.在［a，b］上均匀分布的抽样

在［a，b］上均匀分布的分布函数为：

则例7.β分布

β分布为连续型分布，作为它的一个特例是：其分布函数为：

则例8.指数分布指数分布为连续型分布，其一般形式如下：其分布函数为：

则因为1－ξ也是随机数，可将上式简化为连续性分布函数的直接抽样方法对于分布函数的反函数存在且容易实现的情况，使用起来是很方便的。但是对于以下几种情况，直接抽样法是不合适的。分布函数无法用解析形式给出，因而其反函数也无法给出。分布函数可以给出其解析形式，但是反函数给不出来。分布函数即使能够给出反函数，但运算量很大。下面叙述的挑选抽样方法是克服这些困难的比较好的方法。挑选抽样方法

为了实现从己知分布密度函数f(x)抽样，选取与f(x)取值范围相同的分布密度函数h(x)，如果则挑选抽样方法为：>

即从h(x)中抽样xh，以的概率接受它。下面证明xf

服从分布密度函数f(x)。证明：对于任意x

使用挑选抽样方法时，要注意以下两点：选取h(x)时要使得h(x)容易抽样且M的值要尽量小。因为M小能提高抽样效率。抽样效率是指在挑选抽样方法中进行挑选时被选中的概率。按此定义，该方法的抽样效率E为：所以，M越小，抽样效率越高。

当f(x)在［0，1］上定义时，取h(x)=1，Xh=ξ，

此时挑选抽样方法为>例9.圆内均匀分布抽样令圆半径为R0，点到圆心的距离为r，则r的分布密度函数为分布函数为容易知道，该分布的直接抽样方法是由于开方运算在计算机上很费时间，该方法不是好方法。下面使用挑选抽样方法：取则抽样框图为

＞≤

显然，没有必要舍弃ξ1＞ξ2的情况，此时，只需取就可以了，亦即另一方面，也可证明与具有相同的分布。复合抽样方法

在实际问题中，经常有这样的随机变量，它服从的分布与一个参数有关，而该参数也是一个服从确定分布的随机变量，称这样的随机变量服从复合分布。例如，分布密度函数是一个复合分布。其中Pn≥0，n=1，2，…，且

fn(x)为与参数n有关的分布密度函数，n=1，2，…，参数n服从如下分布

复合分布的一般形式为：其中f2(x/y)表示与参数y有关的条件分布密度函数，F1(y)表示分布函数。 复合分布的抽样方法为:首先由分布函数F1(y)或分布密度函数f1(y)中抽样YF1或Yf1，然后再由分布密度函数f2(x/YF1)中抽样确定Xf2(x/YF)

证明：所以，Xf所服从的分布为f

(x)。例10.指数函数分布的抽样

指数函数分布的一般形式为：引入如下两个分布密度函数：

则使用复合抽样方法，首先从f1(y)中抽取y

再由f2(x/YF1)中抽取x

复合挑选抽样方法

考虑另一种形式的复合分布如下：其中0≤H(x,y)≤M，f2(x/y)表示与参数y有关的条件分布密度函数，F1(y)表示分布函数。抽样方法如下：>

证明：抽样效率为：E=1/M为了实现某个复杂的随机变量y的抽样，将其表示成若干个简单的随机变量x1，x2，…，xn

的函数得到x1，x2，…，xn

的抽样后，即可确定y的抽样，这种方法叫作替换法抽样。即替换抽样方法例11.散射方位角余弦分布的抽样散射方位角φ在［0,2π］上均匀分布，则其正弦和余弦sinφ和cosφ服从如下分布：

直接抽样方法为：

令φ=2θ，则θ在［0,π］上均匀分布，作变换 其中0≤ρ≤1，0≤ρ≤π，则

(x,y)表示上半个单位圆内的点。如果(x,y)在上半个单位圆内均匀分布，则θ在［0,π］上均匀分布，由于因此抽样sinφ和cosφ的问题就变成在上半个单位圆内均匀抽样(x,y)的问题。为获得上半个单位圆内 的均匀点，采用挑选法，在 上半个单位圆的外切矩形内 均匀投点（如图）。 舍弃圆外的点，余下的就是所要求的点。 抽样方法为： 抽样效率

E=π/4≈0.785>为实现散射方位角余弦分布抽样，最重要的是在上半个单位圆内产生均匀分布点。下面这种方法，首先在单位圆的半个外切正六边形内产生均匀分布点，如图所示。于是便有了抽样效率更高的抽样方法： 抽样效率>≤例12.正态分布的抽样标准正态分布密度函数为：引入一个与标准正态随机变量X独立同分布的随机变量Y，则（X,Y）的联合分布密度为：作变换

则（ρ,φ）的联合分布密度函数为： 由此可知，ρ与φ相互独立，其分布密度函数分别为 分别抽取ρ，φ

：

从而得到一对服从标准正态分布的随机变量X和Y：

对于一般的正态分布密度函数N(μ,σ2)的抽样，其抽样结果为：例13.β分布的抽样β分布密度函数的一般形式为：其中n，k为整数。为了实现β分布的抽样，将其看作一组简单的相互独立随机变量的函数，通过这些简单随机变量的抽样，实现β分布的抽样。设x1，x2，…，xn

为一组相互独立、具有相同分布F(x)的随机变量，ζk为x1，x2，…，xn

按大小顺序排列后的第k个，记为：

则ζk的分布函数为：

当F(x)=x

时，

不难验证，ζk的分布密度函数为β分布。因此，β分布的抽样可用如下方法实现： 选取n个随机数，按大小顺序排列后取第k个，即随机抽样的一般方法

加抽样方法

减抽样方法乘抽样方法乘加抽样方法乘减抽样方法对称抽样方法积分抽样方法加抽样方法

加抽样方法是对如下加分布给出的一种抽样方法：其中Pn≥0，

，且

fn(x)为与参数n有关的分布密度函数，n=1，2，…。由复合分布抽样方法可知，加分布的抽样方法为:首先抽样确定n’，然后由fn’(x)中抽样x，即：例14.多项式分布抽样

多项式分布密度函数的一般形式为：将f(x)改写成如下形式：则该分布的抽样方法为：例15.球壳内均匀分布抽样

设球壳内半径为R0，外半径为R1，点到球心的距离为r，则r的分布密度函数为分布函数为 该分布的直接抽样方法是

为避免开立方根运算，作变换： 则x∈[0,1]，其分布密度函数为： 其中

则x及r的抽样方法为：≤≤>>减抽样方法

减抽样方法是对如下形式的分布密度所给出的一种抽样方法：其中A1、A2为非负实数，f1(x)

、f2(x)均为分布密度函数。减抽样方法分为以下两种形式：以上两种形式的抽样方法，究竟选择哪种好，要看f1(x)

、f2(x)哪一个容易抽样，如相差不多，选用第一种方法抽样效率高。

（1）将f

(x)表示为令m表示f2(x)／f1(x)的下界，使用挑选法，从f1(x)中抽取Xf1

抽样效率为：>

（2）将f

(x)表示为使用挑选法，从f2(x)中抽取Xf2

抽样效率为：>例16.β分布抽样

β分布的一个特例： 取A1＝2，A2＝1，f1(x)＝1，f2(x)＝2x，此时m＝0，则根据第一种形式的减抽样方法，有 或＞≤＞≤

由于1－ξ1可用ξ1代替，该抽样方法可简化为： 对于ξ2＞ξ1的情况，可取Xf＝ξ1

，因此 与β分布的推论相同。＞≤

如下形式的分布称为乘分布：其中H(x)为非负函数，

f1(x)为任意分布密度函数。 令M为H(x)的上界，乘抽样方法如下：抽样效率为：乘抽样方法≤＞例17.倒数分布抽样

倒数分布密度函数为： 其直接抽样方法为： 下面采用乘抽样方法，考虑如下分布族： 其中i=1，2，…，该分布的直接抽样方法为：

利用这一分布族，将倒数分布f(x)表示成:

其中， 乘法分布的抽样方法如下:

该分布的抽样效率为：＞≤例18.麦克斯韦(Maxwell)分布抽样

麦克斯韦分布密度函数的一般形式为： 使用乘抽样方法，令 该分布的直接抽样方法为：

此时 则麦克斯韦分布的抽样方法为:

该分布的抽样效率为：＞≤

在实际问题中，经常会遇到如下形式的分布：其中Hn(x)为非负函数，fn(x)为任意分布密度函数，n=1，2，…。不失一般性，只考虑n=2的情况：

将

f(x)改写成如下的加分布形式：乘加抽样方法其中

乘加抽样方法为：

该方法的抽样效率为：>>>≤

这种方法需要知道P1的值(P2=1－P1），这对有些分布是很困难的。下面的方法可以不用计算P1

：

对于任意小于1的正数P1

，令P2=1－P1

；则采用复合挑选抽样方法，有：

当取时，抽样效率最高这时，乘加抽样方法为：>>>≤

由于可知第一种方法比第二种方法的抽样效率高。例19.光子散射后能量分布的抽样令光子散射前后的能量分别为

和（以m0c2

为单位，m0为电子静止质量，c

为光速），，则x

的分布密度函数为： 该分布即为光子散射能量分布，它是由著名的Klin－Nishina

公式确定的。其中K(α)为归一因子：

把光子散射能量分布改写成如下形式： 在［1,1+2α］上定义如下函数:

则有 使用乘加抽样方法：

光子散射能量分布的抽样方法为：

该方法的抽样效率为：＞＞＞≤≤≤

乘减分布的形式为： 其中H1(x)、H2(x)为非负函数，f1(x)、f2(x)为任意分布密度函数。与减抽样方法类似，乘减分布的抽样方法也分为两种。乘减抽样方法

（1）将f

(x)表示为 令H1(x)的上界为M1，的下界为m，使用乘抽 样方法得到如下乘减抽样方法：>

（2）将f

(x)表示为 令H2(x)的上界为M2，使用乘抽样方法，得到另一种乘减抽样方法：>

对称分布的一般形式为：其中f1(x)为任意分布密度函数，满足偶函数对称条件，H(x)为任意奇函数，即对任意x满足：对称分布的抽样方法如下：取η=2ξ－1对称抽样方法>≤

证明： 因为η=2ξ－1，η≤x

相当于ξ≤，因此如下形式的分布密度函数 称为积分分布密度函数，其中f0(x，y)为任意二维分布密度函数，H(x)为任意函数。该分布密度函数的抽样方法为：积分抽样方法>

证明：对于任意x

例22.各向同性散射方向的抽样为了确定各向同性散射方向，根据公式：对于各向同性散射，cosθ在［－1,1］上均匀分布，φ在［0,2π］上均匀分布。由于 直接抽样需要计算三角函数和开方。

定义两个随机变量：

可以证明，当时，随机变量x

和y

服从如下分布:

定义区域为：

则w＝cosθ

的分布可以用上述分布表示成积分分布的形式：

令，则属于上述积分限内的y

一定满足 条件。

各向同性散射方向的抽样方法为： 抽样效率为：＞≤随机抽样的其它方法

偏倚抽样方法近似抽样方法近似-修正抽样方法多维分布抽样方法指数分布的抽样第四节MC方法在材料科学中的应用研究思想发展现状MC方法的具体模拟步骤MC方法模拟晶粒生长的基本原理晶粒长大的基本理论HAZ晶粒长大模拟1.研究思想

蒙特卡罗方法适用于研究材料中的随机过程及现象，该方法将晶体材料分为小的体积单元，将每个体积单元当作一个小单晶，与取向数字相对应，在两个小单晶间赋予能量，通过使系统总能量最小化来完成结构演化模拟。2.发展现状及存在问题

蒙特卡罗方法适用于研究材料中的随机过程及现象。主要应用于模拟薄膜生长、晶粒长大、扩散、相变、缺陷行为、碰撞和渗流等过程。1.以Srolovita，Andeson，Grest和Sahni为主的EXON组最早利用MC方法模拟单相多晶体的晶长大和合金中有序-无序畴的演化;

现在已将此方法推及至有第二相粒子的晶粒长大、异常晶粒长大、二元合金组织演变、静态、动态再结晶等模拟过程。2.发展现状及存在问题

2.国内学者

1994年，陈礼清等利用蒙特卡罗方法模拟了二维多晶体晶粒生长规律。

1995年，宋晓艳等利用三维的蒙特卡罗技术模拟了较完整的单相材料正常晶粒生长的过程，获得了晶粒生长动力学和拓扑学的全面信息，逼真再现了晶粒生长过程。

1999年，钟晓征等以蒙特卡罗方法为基础，对正常晶粒生长形貌演化进行了可视化研究2.发展现状及存在问题

采用Monte-Carlo方法存在以下几方面的问题

（1）点阵离散化导致模拟的晶粒长大形态呈现各向异性；（2）难以直接跟踪生长过程中动态界面的演变；（3）模拟过程只计及最近邻格点的相互作用，难于考虑畸变场等物理场以及晶界形态对晶粒长大过程的影响，不利于在唯象层次

上深入开展对晶粒长大过程的研究。3.MC方法的具体模拟步骤

MC方法的具体步骤如下:(l)在相空间中规定一个初始点--类似于分子动力学和布朗运动中的初始化。(2)产生一个新状态——随机的选取一个新状态或新位形(3)计算跃迁的概率(4)产生一个均匀随机数R(在0与1之间)(5)如果跃迁的概率W小于随机数R，那么把原状态算作一个新状态并回到第2步(6)否则接受新状态并回到第2步

做出判断，是接受还是拒绝一次MC变动4.模拟晶粒长大的基本思想

采用Monte-Carlo模拟晶粒长大的基本思想（1）首先对选取的模拟区域进行离散化。根据模拟区域的不同，将模拟区域离散为二维或三维的微元，每个微元的中心作为离散网格的一个格点。在二维离散化时通常采用三角形、四边形或六边形微元(如图1所示)；在三维离散化时通常采用简单立方体微元；（2）在模拟区域离散化之后，为离散网格的每一个格点随机分配一个从1到Q

之间的晶粒取向数(Q为晶粒取向总数)。

晶粒边界定义为位于具有不同晶粒取向数的两个相邻格点之间，当两个相邻格点的晶粒取向数相同时，则认为此相邻格点为同一个晶粒的一部分。如图2、3所示。图3二维离散图中数字为晶粒取向数，实线为晶粒边界；晶粒边界能量由近邻格点的相互作用指定，近邻格点的相互作用能由Hamilton

方程计算得到，Hamilton方程为：

4.模拟晶粒长大的基本思想

Hamilton方程为：式中：J—设定晶粒边界能量的常数；

Si

—随机选择的格点i的晶粒取向数；

Sj

—格点i的近邻格点j的晶粒取向数；

N—近邻格点总数；

δ—Kronecker的delta函数，当Si与Sj的值相同时，δ取值为1;

当Si与Sj的值不同时，δ取值为0；

4.模拟晶粒长大的基本思想

(3)模拟晶粒边界迁移

MC法晶粒生长过程模拟通过模拟晶粒边界迁移实现。

MC法晶粒边界迁移模拟的实现：

1）随机选择一个格点并将其晶粒取向数重取向为近邻格点的晶粒取向数，由于格点的重取向导致近邻格点的相互作用能改变，重取向以一定的概率被接受。重取向概率定义为：式中ΔE—重取向后近邻格点的相互作用能改变值；

kB—为玻尔兹曼常数；

T—绝对温度；

4.模拟晶粒长大的基本思想

2）格点重取向成功时，晶粒边界发生移动。相反，晶粒边界不发生移动。随机选择的格点总数为离散网格的格点总数时，为一个MC模拟时间步(MCS)。5.晶粒长大的基本理论

1.晶粒长大动力学

2.正常晶粒长大

3.异常晶粒长大（也称为重结晶）晶粒长大动力学

晶粒长大过程中，小晶粒减少而大晶粒继续生长，其结果是系统总的晶界面积逐渐减小，总的界面自由能逐渐降低。所以，界面自由能的降低为晶粒长大提供了驱动力，假设晶粒边界的平均曲率半径与晶粒平均直径d成比例，则晶粒生长的近似平均驱动力为:Ashby和Easterling给出的平均晶粒长大速度为:式中，γgb为晶界gb的表面能，M为晶界迁移率，是温度的函数。由上式可知，晶粒长大驱动力和晶粒长大速度与晶粒直径成反比，晶粒直径越小，晶粒长大的速度越快。如果晶粒长大是在一个极高的温度且没有其它因素抑制晶界的迁移，则一个多晶体材料最终将发展成为一个单晶体。正常晶粒长大正常晶粒长大1、特点：晶粒尺寸一致增长，归一化晶粒尺寸F(R/)和拓扑分布函数P（Nc）不随时间改变其中R为晶粒半径，为平均晶粒半径，Nc是晶粒边数2、生长规律或式中，t为时间，K、B为常数。N≤0.5。注意：当t较大时，若R》R0(初始晶粒半径)时，m=1/n

晶粒的初始结构-Voronoi网格Voronoi网格：

金属材料的每个晶粒都是由众多原子组成的颗粒，这些颗粒的边界是原子在不同的晶向之间排列而形成的。在进行计算机模拟时，这些在晶粒内和晶界上的原子均可用带有晶向的“格点”代替。这种用“格点”表示初始晶粒而形成的网格，即为Voronoi

网格如图所示。

图Voronoi

网格

晶粒的初始结构-Voronoi网格Voronoi网格产生方法：

晶粒生长法、几何法等晶粒生长法具体过程：

1）在平面上随机选取n个格点作为n个晶粒的晶核；这些点可用机内的随机发生器产生，也可用自己编写的程序产生。每个点的坐标为(xi，yi)，2n个随机数便形成了n个点。在晶粒生长法小，这些点便是n个晶核。

2）将n个晶核分别赋以不同的晶向；

3）其余格点均赋以晶向-1，表示未定向；

图Voronoi

网格。

晶粒的初始结构-Voronoi网格4）以n个晶核为核心，每个晶核依次以相同的速度向周围生长（把与该晶粒相邻的格点赋以该晶粒的晶向值，并入该晶粒，使格点成为该晶粒的一部分），使各自的面积越来越大。当这些晶粒的面积增大到互相接触时便停止生长，相互接触处则自动构成各晶粒的边界，直到整个平面被各晶粒占满为止，便构成了n个晶粒的Voronoi

网格。模拟区域中的晶粒数n与真实区域中的晶粒数应相等。

Voronoi

网格特性：一个晶粒内的所有点比其它晶粒内的任意一点都更靠近自己的晶核。晶粒的初始结构-Voronoi网格几何法：一般应用中，几何法比一般晶粒生长法更加常用。产生过程：在平面上随机地撒下n个点后，作每一对点的连线的垂直平分线，这此垂直平分线的交点即为网格中多边形晶粒的顶点，只有那些最靠近晶核的线段构成围绕该核的多边形。图Voronoi

网格算法

晶粒的拓扑结构现象：在晶粒生长的过程中，不可避免地出现大晶粒的长大以牺牲小晶粒为代价(通常所说的“吞噬)的现象，这样，程序中将涉及到一些拓扑变化。晶粒结构的演化包含的基本的拓扑类型:如图所示

(a)双相邻转变变化．即两个晶粒获得一条边，而两个晶粒减少一条边。

(b)三边原胞的消失，即三个相邻原胞减少一条边。其它的变化均可由上述两种变化的组合来实现。晶粒的拓扑结构异常晶粒长大-动力学异常晶粒长大1、特点：在重结晶的显微结构中，一些晶粒的尺寸迅速增大，最大尺寸的晶粒以比算数平均速率大的多的速率增长。2、生长规律

其中X为二次取向的晶粒的面积分数，g(t)为与时间有关的函数。

通常g(t)=apn

其中b,p为常数，大部分p的值为1.8±0.3异常晶粒生长模拟理论基础

在正常晶粒生长的模型基础上调整正常晶粒生长的模拟过程：

根据D.J.Sorolovitz,M.P.Anderson等的做法，将显微结构绘制成分立的晶格(可采用三角、立方或六角晶格)每个格点分配一个1—Q的数字，该数字与所嵌入的晶粒的取向对应，晶界定义为具有不同取向的两个格点的公共边。晶界能由最近邻格点的相互作用来表示。定义晶粒间的这种相互作用为运用Metropolis算法来实现晶粒生长的模拟。异常晶粒生长模拟异常晶粒生长的模拟原理：

和正常晶粒生长的系统，但是体系能量的表达式不同。（原因：体系能量的改变是引起异常晶粒生长的最直接原因，而这种改变与许多因素紧密相关。诸如杂质、应力、对材料进行择优取向等因素都将使生长环境变复杂，从而带来生长的各向异性）。

当材料结构不是很复杂时，通常认为边界迁移率和驱动力对边界运动速度的影响主要，可表示为

V=MP

其中V是边界运动速度，M是边界迁移率，P是驱动力。

通常认为晶粒的异常生长主要是由晶界迁移率各向异性和驱动力各向异性所致，而驱动力，主要考虑晶界能。

重点：由于晶界能和迁移率的各向异性引起的体系能量改变。

异常晶粒生长模拟-物理模型异常晶粒生长的模拟物理模型：1、各向异性的晶界能带来的影响三种处理方式（1）G.S.Gret和D.J.Srolovitz的处理方法体系能量其中θij=2π(sisj)/Q

异常晶粒生长模拟-物理模型（2）A.D.Rollet利D.J.Srolovitz的处理方法按取向将晶粒分为两类：

1）Si≤C，认为晶粒属于类型1；

2）Si〉C，认为晶粒属于类型2。

Si为格点i上Q个取向的一个，

C为区分类型1和类型2而设立的一个取向常数。

晶粒间晶界能为两种（情况A和情况B）

1）情况A

相同晶粒间的能量高(1—1或2—2边界)，不同晶粒间的能量低(1—2边界)，如图所示

异常晶粒生长模拟-物理模型相互作用能其中J1和J2为正常数，且J1>J2。异常晶粒生长模拟-物理模型2）情况B

认为只有2-2间具有较高的能量，如图所示，其相互作用能为：

可以设置不同的C、Q、J1/J2等值来观察显微结构的变化。异常晶粒生长模拟-物理模型（3）Y.Satio处理方法认为晶界能分为两类：高能组（对应大角度晶界）；低能组（对应小角度晶界）引入代表晶界能各向异性的参数γ（高能边界和低能边界的能量比）

0≤γ≤1（