集團變分法(CVM)計算合金的熱力學(xué)量

《計算材料學(xué)》課程設(shè)計指導(dǎo)老師：江建軍 教授電子科學(xué)與技術(shù)系2004年6月1集團變分法（CVM）計算合金的熱力學(xué)量薛旭輝馮鵠顏鑠清馬文武張姣陽胡錦華馬暉葉群余戴君楊昱徐晶摘要：介紹了集團變分法并利用集團變分法建立模型,計算了二元 bcc結(jié)構(gòu)無序合金的巨勢及相對化學(xué)勢等熱力學(xué)量，并用計算機模擬了巨勢及相對化學(xué)勢的變化，摩爾自由能和摩爾熵的變化，得出了比較正確的曲線。關(guān)鍵詞：集團變分法概率變量巨勢相對化學(xué)勢熵自然迭代法TheapplicationofclustervariationmethodincalculationofthermodynamicsquantitiesofalloyAbstract:TheClusterVariationMethod(CVM)whichthemodelaccordingtohavebeenrecommended.Thegrandpotentialoppositechemicalpotentialandsomeotherthermodynamicsquantitiesofdualisticdisorderalloywiththestructureofthebcclatticehavebeencalculated.Thegraphofthevariationofgrandpotentialaccompaniedwiththechangeoftheoppositechemicalpotentialandthevariationofmolefreeenergyaccompaniedwiththechangeofmoleentropyhavebeenchartedbythecomputer.Keywords: cluster variation method; probability variation; grand potential;oppositechemicalpotential;entropy;naturaliterationmethod實驗證明有許多材料,如各種二元合金，金屬間化合物 [1]，插層化合物等在一定條件下發(fā)生有序-無序轉(zhuǎn)變，并且這種轉(zhuǎn)變對材料的理化性能有直接影響。因此，從理論上研究有序-無序相變規(guī)律和相變機制一直是人們關(guān)注的焦點之一。適用于這一要求，Kikuchi 提出了原子集團變分法：（ClusterVariationMethod-CVM），并成功的應(yīng)用于二元合金的有序-無序相變研究。隨后，這一方法及其改進方法被廣泛應(yīng)用于各種有序-無序相變研究，取得了令人矚目的成果。一 原子集團變分法概述有序-無序相變的理論研究始于計算合金的有序-無序相變圖，當(dāng)時均采用Bragg-Williams 近似，即為原子集團變分法 [2]中的點近似。除了一些特殊的情況外，點近似2所推導(dǎo)出的結(jié)果不準(zhǔn)確，只能作一些定性分析。后來，Bethe 考慮了原子之間的最近鄰相互作用（theNearest-NeighborInteraction ）,提出了 Bethe近似（即對近似），類似的還有準(zhǔn)化學(xué)近似，對近似在計算 bcc結(jié)構(gòu)的合金系統(tǒng)上取得了一定的成果，但是，當(dāng)將其應(yīng)用到FCC結(jié)構(gòu)的合金系統(tǒng)時，所給出的結(jié)果與實驗結(jié)果偏差較大，并且不能給出任何穩(wěn)定相。直到后來人們提出了四面體近似，才將此問題解決。原子集團變分法[3]最早由Kikuchi 提出，這一方法的關(guān)鍵在于計算機系統(tǒng)構(gòu)型熵的近似表達式。通過集團組合的階乘推導(dǎo)出了一系列構(gòu)型熵的近似表達式，例如：當(dāng)選擇最近鄰對作為基本原子集團時，系統(tǒng)的構(gòu)型熵為：3 2 2sm=kNA[(2ω-1)∑xilnxi-ω∑∑yij]（1）i=1i=1j=1上式中，k為玻耳茲曼常量，2ω為配位數(shù)，NA為系統(tǒng)中的格點數(shù)原子集團變分法解決相變問題的一般過程如下：1選擇和確定基本原子集團。通過將不同種類的原子放在基本原子集團的各個格點上，來構(gòu)建所有基本原子集團構(gòu)型和次原子集團構(gòu)型。對每一個原子集團構(gòu)型賦予一個組態(tài)變量。根據(jù)這些組態(tài)變量近似寫出整個系統(tǒng)的自由能（或巨勢）的表達式，通常自由能（或巨勢）的表達式將決定近似的精確程度。讓自由能（或巨勢）對組態(tài)變量取極小值，得到一組聯(lián)立代數(shù)方程組，通常方程組是高階代數(shù)方程組。解聯(lián)立代數(shù)方程組，由此可確定系統(tǒng)的平衡態(tài)和相應(yīng)的物理量。在以上五步中，第４步和第5步是最為關(guān)鍵的。一般說來，原子集團選取越大，所得的結(jié)果就越可靠，但是計算量也將呈指數(shù)增長，所以只要選取的原子集團能夠正確反映系統(tǒng)的構(gòu)型就可以了。以前，由于第5步（也就是聯(lián)立代數(shù)方程組）的困難，使這一方法的廣泛應(yīng)用受到限制，但是隨著高性能電子計算機的飛速發(fā)展，這一問題正逐漸被解決。數(shù)學(xué)中求解聯(lián)立代數(shù)方程組，通常有兩種方法：一種是Newton-Raphson迭代法。這種方法需要微分計算和矩陣反演，對迭代初始值的品質(zhì)較敏感，但是收斂較快。另一種是自然迭代法（NaturalIterationMethod），是由Kikuchi提出的。在這個方法中極小值是通過考慮最大集團的占據(jù)概率來獲得的，并且把相關(guān)函數(shù)看作是獨立變量。這種方法避開了微分計算和矩陣反演，對3人為誤差不敏感，只要輸入的起始值能夠反映系統(tǒng)的對稱性，就可以有一個離平衡值很遠的初始值開始迭代，并且能夠收斂，但收斂較慢，計算量較大。二 數(shù)學(xué)模型我們研究一種BCC結(jié)構(gòu)的二元合金[4,5]。在常溫狀態(tài)下，計算其熱力學(xué)量 [6]。為了簡化計算，我們假設(shè)該二元合金無序替位式無限固溶體。圖1 合金結(jié)構(gòu)示意圖圖1中BCC無序相的最近鄰原子對一端在α亞點陣上，另一端在β亞點陣上。表1二元合金無序相最近鄰所需概率變量配置X表示在α格點上出現(xiàn)i原子的概率Y 表示亞點陣上為 i原子而β亞點陣上為j原子的“原子對”概率4點歸一化條件分別為：2∑Xiα=1（2）i=12∑Xiβ=1（3）i=1點概率與對概率的關(guān)系2i=1,2（4）Xiα=∑yijαβj=12j=1,2（5）Xiβ=∑yαβjik=1“對概率”歸一化條件22（6）∑∑y12αβ=1i=1j=1無序相滿足y12=y21 （7）三 熱力學(xué)量的計算摩爾自由能確定系統(tǒng)是否處于平衡狀態(tài)的熱力學(xué)函數(shù)是自由能。CVM模型在定義了集團變量之后，利用這些變量所描述的是Helmholtz自由能。由N（Avogadro常數(shù)）個晶格結(jié)點所構(gòu)成的系統(tǒng)的摩爾Helmholtz自由能為F，熱力學(xué)能（內(nèi)能）為U，熵為S因此有Fm=Um-TSm（8）2熱力學(xué)能（內(nèi)能）關(guān)于熱力學(xué)能，最廣泛接受的表示方法是只計算最近鄰原子對的結(jié)合能[7]的總和，這里也做加處理，即22（9）ωμyijUm=NA∑∑ijij5式中ω是配位數(shù),μij是i和j原子分別處于最近鄰“對”兩端時的最近鄰原子對結(jié)合能。所以應(yīng)該有μij=μji（10）將由上述方法求出的熱力學(xué)能Um和配制熵Sm代入式（8），可以得到Helmholtz自由能Fm與Bragg-Williams近似不同，在集團變分法中不再有通用的Fm表達式，該表達式會因所選擇的原子集團的概率變量而異?！皩啤钡哪朒elmholtz自由能Fm的表達式為22?222?（11）（）lnXi-ω∑∑yijlnyij?Fm=ωNA∑∑μijyij-TkNA?2ω-1∑Xii=1j=1?i=1i=1j=1?求出使這一函數(shù)Fm為極小值的概率變量數(shù)值的過程，也就是求單相平衡狀態(tài)的過程，這個過程叫做變分。由這樣的概率變量數(shù)值所描述的狀態(tài)，才是平衡狀態(tài)。巨勢與相對化學(xué)勢的定義化學(xué)勢是人們熟知的概念，它與摩爾Gibbs自由能的關(guān)系可以容易得到，因而能利用摩'爾自由能曲線加以圖示。相對化學(xué)勢( μi)與化學(xué)勢（μi）有所不同，定義 n元系各組元的相對化學(xué)勢之間的關(guān)系為2'（12）μ=0∑ii=1n元系巨勢GP是由Helmholtz'之間的關(guān)系定義的自由能及F相對化學(xué)勢μiG=F-n''（13）nμmiiPi=1式中ni是組元i的物質(zhì)的量，將式（13）兩邊都除以總的物質(zhì)的量∑可以得到摩爾巨勢GP：ijmnGpm=Fm-∑xiu'i（1i=1（Xi是組元i的摩爾分?jǐn)?shù)）n元系摩爾分?jǐn)?shù)的歸一化條件為：6n∑Xi=1（15）i=1變分與數(shù)值計算變分（Variation）就是求出概率變量數(shù)值的過程，而這個數(shù)值能使熱力學(xué)函數(shù)變成極小值，在保持“相對化學(xué)勢”μi為一定值的前提下，能夠使巨勢Gpm為極小值的狀態(tài)就是平衡態(tài)。所以變分的過程也就是求平衡態(tài)的過程。摩爾巨勢可以表示成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計算合金的熱力學(xué)量的核心思想：使相對化學(xué)勢u'i（注意相對化學(xué)勢u'i和化學(xué)勢ui不同）取定值進行變分?巨勢Gp取極小值時的概率變量(xi,yij),此時得到點(u'i,Gpm),m改變相對化學(xué)勢ui的取值，即可得到關(guān)系曲線Gpm=f(u'i)（18）說明：xi為任意結(jié)點上出現(xiàn)第 i種原子的概率yij為任意相鄰結(jié)點對上出現(xiàn) i-j “原子對”的概率已知配位數(shù)2ω，阿佛加德羅常數(shù) NA，i-j 原子分別處于最近鄰“對”兩端時的原子對結(jié)合能uij，Boltzmann常數(shù)k，熱力學(xué)溫度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，其結(jié)果以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）可推出關(guān)系式： g=

λ

22）kT自此我們可以由以上簡化關(guān)系式（22）利用計算機進行數(shù)值計算，即求解（22）關(guān)于個非線性方程組。

yij的四我們采用自然迭代法（NIM，NaturalIterationMethod ）進行數(shù)值計算。其優(yōu)點是對初值的依賴性小，我們可以粗略進行初值的選擇，因其具有絕對收斂性。自然迭代法步驟：1、一組u'i2、先估算一個適當(dāng)?shù)某跏贾?，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，重復(fù)步驟2—7。得到曲線（u'i，Gpm）四 模擬結(jié)果及分析

（25）（26）（27）摩爾巨勢—相對化學(xué)勢曲線由前面的分析，可以求出摩爾巨勢與相對化學(xué)勢的關(guān)系，進而可以得到熵與溫度的關(guān)系及自由能的關(guān)系。根據(jù)熱力學(xué)的知識，巨勢與相對化學(xué)勢呈現(xiàn)拋物線關(guān)系。在給定的溫度T=300K下，根據(jù)上述介紹的方法，用計算機模擬可以得到右圖所示的關(guān)系。可見當(dāng)相對化學(xué)勢為0時，摩爾巨勢有最小值。 圖2 摩爾巨勢—相對化學(xué)勢摩爾自由能—成分曲線及摩爾自由能—溫度曲線根據(jù)摩爾巨勢—相對化學(xué)勢關(guān)系，可以求出自由能—成分曲線（T=300K） 。模擬結(jié)果如圖3，可見當(dāng)二元相無序合金成分相等時 Fm有最小值。在相對化學(xué)勢取 48kJ/mol，合金各相各占1/2時的Fm—T曲線如圖4所示，隨溫度上升而下降，這與陶瓷的燒結(jié)是吻9合的，且在常溫范圍內(nèi)下降較快。圖 3 摩 爾 自 由 能 — 成 分 曲 線 (300k) 圖 4 摩 爾 自 由 能 — 溫 度 曲 線摩爾熵—成分曲線及摩爾熵—溫度曲線根據(jù)Fm=Um-T*Sm則,可以根據(jù)模擬過程 2得到摩爾熵與成分及溫度的關(guān)系，分別如圖5、圖6所示。由圖5可以看出T=300K，當(dāng)兩相成分相當(dāng)時摩爾熵最大，即系統(tǒng)的混亂程度最大，這與熱力學(xué)理論是相符的。在相對化學(xué)勢取48kJ/mol，合金各相各占1/2時的Sm—T曲線如圖6所示，在一定范圍內(nèi)，當(dāng)溫度上升時，Sm上升很快。但是當(dāng)溫度上升到一定的程度時，摩爾熵會趨于平緩（如圖7