集团变分法(CVM)计算合金的热力学量

《计算材料学》课程设计指导老师：江建军 教授电子科学与技术系2004年6月1集团变分法（CVM）计算合金的热力学量薛旭辉冯鹄颜铄清马文武张姣阳胡锦华马晖叶群余戴君杨昱徐晶摘要：介绍了集团变分法并利用集团变分法建立模型,计算了二元 bcc结构无序合金的巨势及相对化学势等热力学量，并用计算机模拟了巨势及相对化学势的变化，摩尔自由能和摩尔熵的变化，得出了比较正确的曲线。关键词：集团变分法概率变量巨势相对化学势熵自然迭代法TheapplicationofclustervariationmethodincalculationofthermodynamicsquantitiesofalloyAbstract:TheClusterVariationMethod(CVM)whichthemodelaccordingtohavebeenrecommended.Thegrandpotentialoppositechemicalpotentialandsomeotherthermodynamicsquantitiesofdualisticdisorderalloywiththestructureofthebcclatticehavebeencalculated.Thegraphofthevariationofgrandpotentialaccompaniedwiththechangeoftheoppositechemicalpotentialandthevariationofmolefreeenergyaccompaniedwiththechangeofmoleentropyhavebeenchartedbythecomputer.Keywords: cluster variation method; probability variation; grand potential;oppositechemicalpotential;entropy;naturaliterationmethod实验证明有许多材料,如各种二元合金，金属间化合物 [1]，插层化合物等在一定条件下发生有序-无序转变，并且这种转变对材料的理化性能有直接影响。因此，从理论上研究有序-无序相变规律和相变机制一直是人们关注的焦点之一。适用于这一要求，Kikuchi 提出了原子集团变分法：（ClusterVariationMethod-CVM），并成功的应用于二元合金的有序-无序相变研究。随后，这一方法及其改进方法被广泛应用于各种有序-无序相变研究，取得了令人瞩目的成果。一 原子集团变分法概述有序-无序相变的理论研究始于计算合金的有序-无序相变图，当时均采用Bragg-Williams 近似，即为原子集团变分法 [2]中的点近似。除了一些特殊的情况外，点近似2所推导出的结果不准确，只能作一些定性分析。后来，Bethe 考虑了原子之间的最近邻相互作用（theNearest-NeighborInteraction ）,提出了 Bethe近似（即对近似），类似的还有准化学近似，对近似在计算 bcc结构的合金系统上取得了一定的成果，但是，当将其应用到FCC结构的合金系统时，所给出的结果与实验结果偏差较大，并且不能给出任何稳定相。直到后来人们提出了四面体近似，才将此问题解决。原子集团变分法[3]最早由Kikuchi 提出，这一方法的关键在于计算机系统构型熵的近似表达式。通过集团组合的阶乘推导出了一系列构型熵的近似表达式，例如：当选择最近邻对作为基本原子集团时，系统的构型熵为：3 2 2sm=kNA[(2ω-1)∑xilnxi-ω∑∑yij]（1）i=1i=1j=1上式中，k为玻耳兹曼常量，2ω为配位数，NA为系统中的格点数原子集团变分法解决相变问题的一般过程如下：1选择和确定基本原子集团。通过将不同种类的原子放在基本原子集团的各个格点上，来构建所有基本原子集团构型和次原子集团构型。对每一个原子集团构型赋予一个组态变量。根据这些组态变量近似写出整个系统的自由能（或巨势）的表达式，通常自由能（或巨势）的表达式将决定近似的精确程度。让自由能（或巨势）对组态变量取极小值，得到一组联立代数方程组，通常方程组是高阶代数方程组。解联立代数方程组，由此可确定系统的平衡态和相应的物理量。在以上五步中，第４步和第5步是最为关键的。一般说来，原子集团选取越大，所得的结果就越可靠，但是计算量也将呈指数增长，所以只要选取的原子集团能够正确反映系统的构型就可以了。以前，由于第5步（也就是联立代数方程组）的困难，使这一方法的广泛应用受到限制，但是随着高性能电子计算机的飞速发展，这一问题正逐渐被解决。数学中求解联立代数方程组，通常有两种方法：一种是Newton-Raphson迭代法。这种方法需要微分计算和矩阵反演，对迭代初始值的品质较敏感，但是收敛较快。另一种是自然迭代法（NaturalIterationMethod），是由Kikuchi提出的。在这个方法中极小值是通过考虑最大集团的占据概率来获得的，并且把相关函数看作是独立变量。这种方法避开了微分计算和矩阵反演，对3人为误差不敏感，只要输入的起始值能够反映系统的对称性，就可以有一个离平衡值很远的初始值开始迭代，并且能够收敛，但收敛较慢，计算量较大。二 数学模型我们研究一种BCC结构的二元合金[4,5]。在常温状态下，计算其热力学量 [6]。为了简化计算，我们假设该二元合金无序替位式无限固溶体。图1 合金结构示意图图1中BCC无序相的最近邻原子对一端在α亚点阵上，另一端在β亚点阵上。表1二元合金无序相最近邻所需概率变量配置X表示在α格点上出现i原子的概率Y 表示亚点阵上为 i原子而β亚点阵上为j原子的“原子对”概率4点归一化条件分别为：2∑Xiα=1（2）i=12∑Xiβ=1（3）i=1点概率与对概率的关系2i=1,2（4）Xiα=∑yijαβj=12j=1,2（5）Xiβ=∑yαβjik=1“对概率”归一化条件22（6）∑∑y12αβ=1i=1j=1无序相满足y12=y21 （7）三 热力学量的计算摩尔自由能确定系统是否处于平衡状态的热力学函数是自由能。CVM模型在定义了集团变量之后，利用这些变量所描述的是Helmholtz自由能。由N（Avogadro常数）个晶格结点所构成的系统的摩尔Helmholtz自由能为F，热力学能（内能）为U，熵为S因此有Fm=Um-TSm（8）2热力学能（内能）关于热力学能，最广泛接受的表示方法是只计算最近邻原子对的结合能[7]的总和，这里也做加处理，即22（9）ωμyijUm=NA∑∑ijij5式中ω是配位数,μij是i和j原子分别处于最近邻“对”两端时的最近邻原子对结合能。所以应该有μij=μji（10）将由上述方法求出的热力学能Um和配制熵Sm代入式（8），可以得到Helmholtz自由能Fm与Bragg-Williams近似不同，在集团变分法中不再有通用的Fm表达式，该表达式会因所选择的原子集团的概率变量而异。“对近似”的摩尔Helmholtz自由能Fm的表达式为22?222?（11）（）lnXi-ω∑∑yijlnyij?Fm=ωNA∑∑μijyij-TkNA?2ω-1∑Xii=1j=1?i=1i=1j=1?求出使这一函数Fm为极小值的概率变量数值的过程，也就是求单相平衡状态的过程，这个过程叫做变分。由这样的概率变量数值所描述的状态，才是平衡状态。巨势与相对化学势的定义化学势是人们熟知的概念，它与摩尔Gibbs自由能的关系可以容易得到，因而能利用摩'尔自由能曲线加以图示。相对化学势( μi)与化学势（μi）有所不同，定义 n元系各组元的相对化学势之间的关系为2'（12）μ=0∑ii=1n元系巨势GP是由Helmholtz'之间的关系定义的自由能及F相对化学势μiG=F-n''（13）nμmiiPi=1式中ni是组元i的物质的量，将式（13）两边都除以总的物质的量∑可以得到摩尔巨势GP：ijmnGpm=Fm-∑xiu'i（1i=1（Xi是组元i的摩尔分数）n元系摩尔分数的归一化条件为：6n∑Xi=1（15）i=1变分与数值计算变分（Variation）就是求出概率变量数值的过程，而这个数值能使热力学函数变成极小值，在保持“相对化学势”μi为一定值的前提下，能够使巨势Gpm为极小值的状态就是平衡态。所以变分的过程也就是求平衡态的过程。摩尔巨势可以表示成m≡mn（16）GpF-xμi=1将式（11）带入式（16）可以得到“对近似摩尔巨势”表达式2 2 2pm≡ωNuy-TkN[(2ω-1)xlnx-ωy]-xμ（17）GA∑∑ijijA∑ii∑∑ij∑iiiji=1i=1j=1iCVM计算合金的热力学量的核心思想：使相对化学势u'i（注意相对化学势u'i和化学势ui不同）取定值进行变分?巨势Gp取极小值时的概率变量(xi,yij),此时得到点(u'i,Gpm),m改变相对化学势ui的取值，即可得到关系曲线Gpm=f(u'i)（18）说明：xi为任意结点上出现第 i种原子的概率yij为任意相邻结点对上出现 i-j “原子对”的概率已知配位数2ω，阿佛加德罗常数 NA，i-j 原子分别处于最近邻“对”两端时的原子对结合能uij，Boltzmann常数k，热力学温度T计算用表达式：222u'i+u'j（19）Gpm=ωNA∑∑uijyij-TkNA[(2ω-1)∑xilnxi-ω∑∑yijlnyij]-∑∑2yijiji=1i=1j=1ij22用Lagrange不定乘子λ的形式将“对近似”归一化条件∑∑yij=1引入Gpm中，将摩尔巨i=1j=1势除以NAkT，其结果以g=Gpm记之。TkNA7uij222u'i+u'jy+λ(1-g=ωy-(2ω-1)xlnx+ωylny-y)ijiiijijijjkT2kTkTijii=1i=1j=1ijij（本式中的u'i，u'j等为上式的1NA）取微分?g=0(i,j=1,2)得：?yij2ω-1uij+u'i+u'j)yij=(xixj)ωexp(λ-(i,j=1,2)ωkTkT2ωkT由于λ~22y=1，且g=22?g∑∑?y∑∑iji=1j=1i=1j=1ij

（20）21）可推出关系式： g=

λ

22）kT自此我们可以由以上简化关系式（22）利用计算机进行数值计算，即求解（22）关于个非线性方程组。

yij的四我们采用自然迭代法（NIM，NaturalIterationMethod ）进行数值计算。其优点是对初值的依赖性小，我们可以粗略进行初值的选择，因其具有绝对收敛性。自然迭代法步骤：1、一组u'i2、先估算一个适当的初始值?，xi3、取过渡量η计算ij2ω-1uij+u'i+u'j)η=(xx)ωexp(-（23）ijijkT2ωkT?4、求取λy=exp(?η=1λ)η=exp(λ)（24）ijijijijijij?8?-1?-1λ∑∑ij?∑∑ijexp(λ=ωkTln()=(ηη)ωkTijij5、阶段的?yij??λijijωkT6、用式xi?，求?=∑yijxijxiλ7、将?迭代计算步骤3—6，直至?不再收敛为止。=λkTGpm=NAkTg得到点（u'i，Gpm）8、另取u'i，重复步骤2—7。得到曲线（u'i，Gpm）四 模拟结果及分析

（25）（26）（27）摩尔巨势—相对化学势曲线由前面的分析，可以求出摩尔巨势与相对化学势的关系，进而可以得到熵与温度的关系及自由能的关系。根据热力学的知识，巨势与相对化学势呈现抛物线关系。在给定的温度T=300K下，根据上述介绍的方法，用计算机模拟可以得到右图所示的关系。可见当相对化学势为0时，摩尔巨势有最小值。 图2 摩尔巨势—相对化学势摩尔自由能—成分曲线及摩尔自由能—温度曲线根据摩尔巨势—相对化学势关系，可以求出自由能—成分曲线（T=300K） 。模拟结果如图3，可见当二元相无序合金成分相等时 Fm有最小值。在相对化学势取 48kJ/mol，合金各相各占1/2时的Fm—T曲线如图4所示，随温度上升而下降，这与陶瓷的烧结是吻9合的，且在常温范围内下降较快。图 3 摩 尔 自 由 能 — 成 分 曲 线 (300k) 图 4 摩 尔 自 由 能 — 温 度 曲 线摩尔熵—成分曲线及摩尔熵—温度曲线根据Fm=Um-T*Sm则,可以根据模拟过程 2得到摩尔熵与成分及温度的关系，分别如图5、图6所示。由图5可以看出T=300K，当两相成分相当时摩尔熵最大，即系统的混乱程度最大，这与热力学理论是相符的。在相对化学势取48kJ/mol，合金各相各占1/2时的Sm—T曲线如图6所示，在一定范围内，当温度上升时，Sm上升很快。但是当温度上升到一定的程度时，摩尔熵会趋于平缓（如图7