集團(tuán)變分法(CVM)計(jì)算合金的熱力學(xué)量

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

yij的四我們采用自然迭代法（NIM，NaturalIterationMethod ）進(jìn)行數(shù)值計(jì)算。其優(yōu)點(diǎn)是對(duì)初值的依賴性小，我們可以粗略進(jìn)行初值的選擇，因其具有絕對(duì)收斂性。自然迭代法步驟：1、一組u'i2、先估算一個(gè)適當(dāng)?shù)某跏贾?，xi3、取過(guò)渡量η計(jì)算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、將?迭代計(jì)算步驟3—6，直至?不再收斂為止。=λkTGpm=NAkTg得到點(diǎn)（u'i，Gpm）8、另取u'i，重復(fù)步驟2—7。得到曲線（u'i，Gpm）四 模擬結(jié)果及分析

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