葉小彬，何志海，潘必才,,*

(1.中國科學技術大學 合肥微尺度物質科學國家研究中心，安徽 合肥 230026；2.中國科學技術大學 物理系，安徽 合肥 230026)

聚變堆的安全可靠運行依賴較多因素，其中，面向等離子體的第一壁材料的穩(wěn)定性是關鍵因素之一[1]。當聚變堆在正常運行時，第一壁材料承受著高能粒子的輻照和高熱流密度的沖擊[2-3]。在如此極端的環(huán)境中，面向等離子體的第一壁材料的結構會發(fā)生急劇的變化，從而強烈影響材料的熱力學性質和微結構?，F(xiàn)已采用實驗和理論模擬對第一壁材料的微結構和熱力學性能進行了廣泛的研究[4-7]。

由固體物理可知，當材料的電子結構發(fā)生變化時，材料的結構和性質也會發(fā)生相應改變。因此，在物理上，為深入研究和理解材料的結構與物理性質的關聯(lián)，需密切關注材料的電子結構。在核聚變堆中服役的材料身處極端復雜的環(huán)境，這些環(huán)境因素會強烈調制材料的電子結構。首先，聚變堆在正常運行的情況下，聚變堆中等離子體芯部的溫度約為100 000 000 K。根據(jù)黑體輻射理論中的維恩位移定律，高溫等離子體會源源不斷地向其周圍的第一壁材料鎢(W)輻射以X射線為主體的高能光子。這些高能光子入射到材料，與材料中的原子作用，激發(fā)甚至電離原子中內殼層的部分電子。于是，第一壁材料W中有大量的電子處于激發(fā)態(tài)。其次，從等離子體芯部會溢出大量的高能粒子。這些粒子既有帶電的離子也有電中性的中子。帶電的離子雖受削刮層的控制，釋放出來后的能量被減弱了，但仍有較高的動能，轟擊第一壁材料后，將部分動能傳遞給被碰撞的W原子。這些獲得較高動能的W原子發(fā)生高速位移。在高速位移時，W原子核外的外層電子因受原子核較弱的吸引作用而脫離W原子，因而此時高速位移的W原子其實就是高速位移的W離子。高速位移的離子通過庫侖作用，誘導其周圍的W原子中的部分電子被激發(fā)。處于激發(fā)態(tài)的電子是不穩(wěn)定的，通常會經(jīng)歷一個特征時間后發(fā)生退激發(fā)。退激發(fā)有多種途徑，如電子從高能級躍遷到低能級，釋放出光子；或電子將部分能量傳遞給晶格，從高能態(tài)轉變成低能態(tài)，這是一個無輻射的躍遷過程。無輻射躍遷包含著電子與聲子的耦合和非熱行為。其中，非熱行為是指通過改變電子態(tài)使離子間的受力發(fā)生改變，從而改變材料的結構。在非熱過程中，激發(fā)的電子將部分能量通過量子力學力做功傳遞給了晶格。在W中，無輻射躍遷是最主要的電子退激發(fā)的方式。

已有實驗對聚變堆中的第一壁材料W受輻照后的微結構等物理問題進行了研究[8-10]?，F(xiàn)有的大量實驗數(shù)據(jù)來自于對輻照后的第一壁材料的探測，此時，材料不處于輻照的環(huán)境，所獲得的數(shù)據(jù)包含著輻照后遺留的效果，不能真實反映輻照時對材料產(chǎn)生的實時效應。理論上也有許多的研究工作探究粒子輻照對W的微結構和熱力學性能的影響[11-13]。在大量的關于核材料的理論模擬中，已有的主體的理論研究工作均在經(jīng)驗分子動力學和經(jīng)典蒙特卡羅的構架中開展，缺少對電子激發(fā)效應的考慮。為在理論模擬中考慮電子激發(fā)效應，提出了雙溫度模型的理論模擬方案[14-16]和含時的緊束縛動力學方法[17]，但前者未考慮真實的量子態(tài)，后者采用單s態(tài)的緊束縛動力學研究離子輻照中的電子阻止本領，獲得了較好的結果，但無法研究輻照時處于電子激發(fā)態(tài)下材料的結構演化規(guī)律，也不能計算材料的其他物理性質。

本文發(fā)展W體系的緊束縛勢模型，在該勢模型中，離子間的相互作用包含量子力學力，同時引入描述電子在不同激發(fā)態(tài)分布的電子熵。將緊束縛勢模型與分子動力學結合，對處于不同激發(fā)態(tài)下材料的微結構演化進行研究。

1 理論研究方法

1.1 緊束縛勢模型中的電子態(tài)

對于處于定態(tài)的體系，可通過求解其單粒子能量本征方程(式(1))獲得其電子結構。

Hψk(r)=εkψk(r)

(1)

式中：H為體系的哈密頓算符；k為布里淵區(qū)的波矢；εk為體系的能量本征值；ψk(r)為相應的本征波函數(shù)。當體系的原子軌道波函數(shù)在空間上比較定域時，體系的本征波函數(shù)可在原子軌道波函數(shù)組合成布洛赫和函數(shù)為基的表象中展開：

(2)

式中：Ri為體系中第i個晶格矢量；φα(r-Ri)為體系中位于第i個晶格中位置矢量為r的原子的α軌道波函數(shù)；Ci,α為對應的展開系數(shù)。將式(2)代入式(1)，可獲得體系的哈密頓矩陣元：

(3)

(4)

由式(3)和(4)可知，對一個體系的哈密頓矩陣的構造轉化為對hjβ,iα的求解。式(3)包含同一個原子內部的庫侖積分項hiα,iα和不同原子間的跳躍積分項。其中，hiα,iα可近似地等于相應的原子軌道能級。跳躍積分項中包含多中心積分，因而hjβ,iα的計算量很大。為簡便計算，采用Slater和Koster提出的雙中心近似[18]替代hjβ,iα中的多中心積分。

在Slater和Koster的雙中心近似下，定義一些不同類型的鍵積分，如對1個只包含s、p原子軌道的體系中2個成鍵的原子i和j，它們之間存在著以下4種類型的鍵積分：

(5)

采用Slater和Koster定義的鍵積分可表達出上述的積分項，具體表達式可參考文獻[18]。

為快速求解大規(guī)模復雜體系，進一步用經(jīng)驗公式表達上述鍵積分，從而避開大量復雜的積分運算。在構造的W體系的緊束縛勢中，采用下列的經(jīng)驗公式表達鍵積分[19]：

(6)

Rij=rij(1+δ1Δ+δ2Δ2+δ3Δ3)

(7)

(8)

(9)

(10)

式(6)中的ll′m為鍵積分的類型，l、l′=s、p或d，m=σ、π或δ。為增強緊束縛勢模型的可傳遞性，還考慮了環(huán)境原子的影響。環(huán)境原子對鍵積分的影響體現(xiàn)在式(9)中環(huán)境原子k對原子i和j之間的有效距離Rij(有效鍵長)上(式(7)，其中rij為原子i與j之間的距離)，環(huán)境原子k的構型決定了Rij的取值。在式(6)～(9)中，α1、α2、α3、α4、δ1、δ2、δ3、β1和β2是需通過擬合過程確定的參數(shù)。式(10)為鍵積分的截斷函數(shù)，即鍵積分函數(shù)在原子之間的有效距離超過rcut后，會快速衰減至0。文獻[19]給出了W體系的緊束縛勢模型的具體參數(shù)。

1.2 緊束縛勢模型中的總能

在緊束縛近似下，體系的總能Etot為：

Etot=Eband+Erep

(11)

(12)

式中：Eband為帶結構能；εl和nl分別為體系中第l能量本征態(tài)對應的本征值和電子占據(jù)數(shù)，εl可通過求解本征值方程得到；fl為Fermi-Dirac分布函數(shù)；Erep為體系的總排斥能，可表達為原子間的排斥相互作用能，寫成原子間對勢的求和：

(13)

為增強緊束縛勢模型的可傳遞性，采用一種多項式表述的方案，即：

(14)

F(x)=C0+C1x+C2x2+C3x3+C4x4

(15)

為減小計算量，緊束縛勢模型通常采用最小基近似。因此在選取原子軌道基組時，只考慮價電子軌道以及最低的空軌道。對于W體系，選取5d、6s以及6p軌道作為原子軌道基組。

此外，還需要構造體現(xiàn)原子間排斥作用的兩體排斥勢，即式(14)中的φ(rij)。為了方便，丁文藝等選取與鍵積分表達式類似的經(jīng)驗公式作為兩體排斥勢的表達式，不同之處在于二者的參數(shù)不同，具體參數(shù)可參考文獻[19]。

本文所使用的W緊束縛勢模型已被用于研究完美W晶體的晶格振動譜、彈性常數(shù)、W中的點缺陷(自間隙原子、單原子空位)形成能和晶格熔化的熔點[19]、W在有限溫度下的電導率和熱導率[20-21]。這些研究結果與相關的實驗結果和第一性原理計算的結果一致。

1.3 電子熵

量子力學的微擾論認為，一個處于定態(tài)的微觀體系受突發(fā)微擾時，體系的定態(tài)不改變。高能光子和高能粒子作用到第一壁材料W時導致材料中的部分電子突然變成了激發(fā)態(tài)電子，而部分激發(fā)態(tài)電子通過非熱效應將部分能量傳遞給晶格也可看成是體系的突發(fā)行為。因此，在研究這類突發(fā)行為時，可使用體系的定態(tài)方程求解，但同時要加入由于部分電子被激發(fā)而引起的電子熵對能量的貢獻。

體系中的電子在能量空間中的分布遵從Fermi-Dirac統(tǒng)計分布。當體系中的一群電子處于不同的激發(fā)態(tài)時，其電子熵為：

(16)

式中，fi為Fermi-Dirac分布函數(shù)。于是，體系的自由能為：

(17)

1.4 緊束縛分子動力學

為在有限溫度下對體系的原子結構的演化進行模擬，將緊束縛勢模型與分子動力學結合。在緊束縛分子動力學(TBMD)中，體系的自由能對原子的位置坐標求偏導數(shù)：

(18)

即可獲得該原子在各坐標方向上所受到的力。將原子所受到的力代入牛頓第二運動定律，采用差分方程求解[22]，可獲得原子在受力時的空間位置隨時間的演化。

對式(18)進行展開：

(19)

式中：rl為原子l的位置矢量；-?Erep/?rl為排斥能貢獻的力；-?Eband/?rl為帶結構能貢獻的力，又稱為Hellmann-Feynman力，其具體的表達式為：

(20)

式中，|ψk〉為電子的能量本征矢量，可表述成體系的原子軌道波函數(shù)的線性組合。

在分子動力學模擬過程中，每一時刻體系全部的N個粒子運動的速度為{vi(i=1,N)}，采用約束條件：

(21)

便能控制體系溫度為Tion(Tion為體系的離子溫度，不是電子溫度Tel，當體系處于平衡態(tài)時，其離子溫度與電子溫度才會相等)時的演化。

由于TBMD中力的表達式包含量子力學力，因此較經(jīng)驗勢，緊束縛勢在描述原子間相互作用時包含經(jīng)驗勢所缺失的內容。

2 結果與討論

2.1 微結構中的電子溫度效應

首先計算在有限溫度的熱力學平衡態(tài)下W晶體的晶格變化。由于體系處于熱力學平衡狀態(tài)，故體系中的晶格溫度與電子溫度相同，所以在計算中，Tion=Tel。圖1a中的藍色曲線示出了當體系的溫度從500 K變化到3 500 K過程中W晶體的晶格膨脹率。計算的晶格變化率隨溫度的變化關系與實驗結果幾乎一致。為揭示晶格溫度和電子溫度在晶格膨脹中的效果，在計算中還分別考慮了Tion=0 K和Tel=0 K的情況。當Tel=0 K時，體系的晶格溫度從500 K逐漸升高到3 500 K，此時，晶格的膨脹率均系統(tǒng)地低于Tion=Tel的情形，如圖1a中綠色曲線所示。而當Tion=0 K時，如圖1b所示，Tel從0 K升高到20 000 K時，體系的晶格常數(shù)也隨Tel的增大而單調增大，且本文中基于緊束縛理論的計算結果與基于密度泛函理論的計算結果十分吻合。特別是當Tel大于5 000 K后，晶格快速膨脹。計算結果表明：在中等溫度(～5 000 K)以下，W晶體的晶格膨脹主要由晶格溫度驅動，但在溫度較高時，電子溫度導致的晶格膨脹效應不可被忽略。特別是當電子溫度很高(>10 000 K)時，即便晶格溫度不高，電子溫度也會導致很大程度的晶格膨脹。電子溫度導致的晶格膨脹是激發(fā)態(tài)電子的非熱效應的一種宏觀表現(xiàn)。激發(fā)態(tài)電子的非熱效應是很抽象的復雜效應，建議采用實時原位測量材料在輻照下的晶格變化來檢驗上述物理效應。

a——體系的溫度500～3 500 K過程中W的晶格膨脹率；b——Tion=0 K下W的晶格常數(shù)隨體系電子溫度的變化圖1 W體系的晶格變化Fig.1 Lattice change in W system

第一壁材料W擁有很大的表面。通常，W表面及其附近區(qū)域的原子受到的勢場與W塊體中原子受到的勢場不同。采用由28個原子層構成的薄層模型模擬W的表面區(qū)域，薄層的表面取向為[001]方向，研究當表面附近區(qū)域承受高能光子輻照時其電子的非熱效果。計算中，選取Tion=1 000 K和Tel=20 000 K進行TBMD模擬，并觀察在模擬過程中一些不同時刻的薄層結構，如圖2所示。

圖2 W體系的結構演化Fig.2 Structural evolution of W system

顯然，隨著電子溫度的加入，薄層迅速膨脹。當TBMD的模擬時間達2.6 ps時，薄層的中間區(qū)域開始呈現(xiàn)無序化。隨著模擬時間的增長，這種無序化的程度加重，并快速轉變成非晶的無序結構。同時，在無序結構的區(qū)域中出現(xiàn)了微孔。這些微孔的形狀和尺度均隨時間演化。這種結構的演化現(xiàn)象與飛秒激光輻照W材料時所呈現(xiàn)的現(xiàn)象非常類似[23]。值得注意的是，在結構演化過程中，靠近表面層的原子結構未呈現(xiàn)強烈的無序化特征，而是幾乎保持原有的晶體結構特征。這與薄層中心區(qū)域的結構演化規(guī)律大相徑庭。

圖3 W的最近鄰原子間相互作用力隨電子溫度的變化Fig.3 Interaction force among the nearest-neighbor atoms varies with electronic temperature

為理解上述的結構演化規(guī)律，分析體系在結構演化過程中原子所受到的力。如前所述，在緊束縛勢模型中，體系中每個原子受到的力由電子結構引起的吸引力和原子間的排斥力構成。圖3示出了W的最近鄰原子間相互作用力的大小隨體系電子溫度的變化關系?？煽吹剑旙w系的電子溫度為20 000 K時，所有原子受到的吸引力均大幅度減小，而電子處于激發(fā)態(tài)時，激發(fā)態(tài)電子對原子間的排斥勢無直接影響。雖然如此，由于原子間的吸引力減小了，原子間的間距就會增大，這也相應地減小了原子間的排斥力。于是，從整體上看，當體系處于電子的高激發(fā)態(tài)時，體系的晶格大幅度膨脹。這種晶格膨脹是由激發(fā)態(tài)電子驅動的非熱效應所導致。

隨著分子動力學模擬時間的增長，薄膜上下表面原子層的原子通過馳豫緩慢釋放應力。當原子層接近薄層的中心區(qū)域，這些原子層上非熱效應導致的力必須通過與鄰近原子層的作用釋放，薄膜中間區(qū)域中的原子受到的非熱力不易快速釋放。再加上在晶格溫度驅動下原子的無規(guī)運動，使得薄膜中間區(qū)域的原子表現(xiàn)出從有序結構向非晶無序結構的演化，并且中間部位的原子受到的非熱應力通過其周圍原子的有效馳豫減小，進而出現(xiàn)了孔洞。模擬時間足夠長后，應力均被釋放，部分孔洞逐漸被復合。

2.2 W表面的預熔化

通常，固體材料的表面熔化行為與塊體熔化的行為不同。在表面的原子層，由于其中原子所處勢場的特殊性，這類原子層熔化所對應的熔點要低于塊體中的熔點，即表面預熔化現(xiàn)象。為研究W表面的預熔化，選擇圖4a所示的具有(110)表面結構的模型。在這個結構模型中，原子層的厚度約為3 nm，它由3個部分構成：底部按體相結構被固定的3個原子層，中間作為緩沖區(qū)的3個原子層和上面用于模擬表面的8個原子層。TBMD模擬中的控溫條件均施加到用于模擬表面的8個原子層的區(qū)域。在各溫度下各雙層的均方鍵長漲落δ為：

(22)

丁文藝等[19]在只考慮晶格溫度效應而忽略電子溫度效應時研究W的表面融化行為，發(fā)現(xiàn)原子層越靠近表面，局部的熔點變得越低。換言之，第一壁材料W的表面會表現(xiàn)出低于塊體的熔點。在本文中，考慮電子溫度效應，即Tion=Tel。其中，溫度從1 000 K逐漸增大到3 500 K，如圖4b所示。顯然，考慮了電子溫度后，W表面的熔點相對于未考慮電子溫度時的熔點有所降低。特別是表面上的雙原子層的熔點僅約3 000 K，低于未考慮電子溫度時的3 500 K[19]。因此，當體系處于電子激發(fā)態(tài)時，表面的熔化溫度也會下降。這種表面預熔化現(xiàn)象對于評估在聚變堆中服役的第一壁材料W的熱力學性能十分重要。

a——分子動力學設置；b——Tion=Tel時不同深度的表面雙原子層的均方鍵長漲落隨溫度的變化圖4 W表面的分子動力學模擬Fig.4 Molecular dynamics simulation of W surface

3 結論

本文介紹了W體系的緊束縛勢模型。將緊束縛勢模型與分子動力學結合，研究W材料的微結構演化與熱力學性質的電子溫度效應。研究發(fā)現(xiàn)第一壁材料W表面的熔點低于塊體的熔點，也低于未考慮電子溫度時W表面的熔點。對W材料在非平衡條件(電子溫度遠高于離子溫度)下的微結構演化進行研究，發(fā)現(xiàn)W的金屬鍵隨體系電子溫度的升高逐漸變弱，導致W晶格的劇烈熱膨脹。此外，當體系瞬間處于足夠高的電子溫度下，劇烈的熱膨脹將導致體系在較低的離子溫度下(遠低于熔點)出現(xiàn)內部熔化而表面區(qū)域不熔化的反?，F(xiàn)象。需強調，由于材料在聚變堆中受到粒子輻照時，被輻照的區(qū)域是局部和隨機的，并非同時均勻發(fā)生在全部的壁材料上。因此，材料中那些被強烈輻照的局部區(qū)域才可能表現(xiàn)出上述預測的現(xiàn)象。

感謝中國科學技術大學超級計算中心提供的計算支持。