實體動力學(xué)蒙特卡羅方法在核材料輻照損傷模擬中的應(yīng)用

潘卓楠，楊章燦,*，周洪波，呂廣宏

(1.華中科技大學(xué)，湖北 武漢 430074；2.北京航空航天大學(xué)，北京 100191)

核反應(yīng)堆中嚴峻的運行條件對反應(yīng)堆裝置材料提出了巨大的挑戰(zhàn)。材料暴露在高輻照水平的系統(tǒng)中，高能粒子的撞擊會破壞材料原本完美的晶格，產(chǎn)生諸如位錯、空洞、析出相等缺陷，而這些缺陷嚴重影響著材料的宏觀性能[1-2]。輻照引起材料內(nèi)部微結(jié)構(gòu)演化在時間和空間上是個多尺度的過程，在實驗上，很難實現(xiàn)原子空間尺度、微納秒時間尺度上對實驗的精準控制來進行缺陷演化機制的研究。而計算機模擬能夠彌補實驗上的不足，可用于解釋微結(jié)構(gòu)演化的機制，并且能夠根據(jù)輻照條件對輻照效應(yīng)提出預(yù)測。

輻照導(dǎo)致微結(jié)構(gòu)演化的多尺度特性很難用單一的計算模型來研究，因此微結(jié)構(gòu)演化模擬的研究通常使用多尺度模擬方法[3]。入射高能粒子與被輻照材料基體原子碰撞后產(chǎn)生初級撞出原子(PKA)，其能量分布可以通過粒子輸運程序，如MCNP程序來計算。隨后，PKA引發(fā)級聯(lián)碰撞(<1 ps)產(chǎn)生大量的點缺陷，這個過程通常導(dǎo)致熱峰的產(chǎn)生，極大升高級聯(lián)區(qū)域的溫度，可認為發(fā)生熔化。熱峰區(qū)的原子不斷將能量傳遞給周圍的原子，由熔化狀態(tài)快速(<10 ps)回到凝聚或淬火冷卻階段，最終建立熱力學(xué)平衡態(tài)。冷卻期間形成了穩(wěn)定的點陣缺陷，包括點缺陷和缺陷團。以上過程通常稱為初級輻照損傷，時間尺度約為10 ps，空間尺度約為數(shù)十納米，可使用二體碰撞近似(BCA)方法進行缺陷數(shù)的快速估算，更為準確的計算通常使用分子動力學(xué)(MD)方法。

初級輻照損傷形成的穩(wěn)定缺陷可發(fā)生熱激發(fā)擴散，擴散到其他區(qū)域與級聯(lián)內(nèi)其他缺陷、材料中原有缺陷或其他級聯(lián)區(qū)中的缺陷發(fā)生互相作用，最終形成空洞、氣泡、位錯等宏觀缺陷。這一過程即微結(jié)構(gòu)演化，時間和空間尺度均已達到實驗觀測的尺度，一般使用的模擬方法有動力學(xué)蒙特卡羅(KMC)、團簇動力學(xué)(CD)和相場方法(PF)。這些模擬方法通常需要MD或密度泛函理論(DFT)計算提供缺陷的熱力學(xué)和動力學(xué)參數(shù)。

本文介紹的KMC方法是模擬受輻照材料微結(jié)構(gòu)演化的有效工具，能夠在保留缺陷空間關(guān)聯(lián)性的前提下處理缺陷的擴散和缺陷之間的相互作用，并且具有足夠大的時空尺度，使其模擬結(jié)果能夠直接與實驗值進行對比。KMC是一種基于馬爾可夫過程的方法，在給定體系中各事件的發(fā)生速率之后，可以用于模擬體系隨時間的演化。在材料模擬中，KMC方法不關(guān)注原子的具體運動軌跡，而將其粗粒化為體系組態(tài)躍遷的過程。組態(tài)的躍遷是無記憶性的，僅與躍遷的概率相關(guān)。KMC方法直接處理狀態(tài)到狀態(tài)的轉(zhuǎn)換，因此相比于MD模擬方法，可以獲得幾個數(shù)量級的時間尺度增益。

根據(jù)對缺陷處理方式的差異，KMC方法可分為原子動力學(xué)蒙特卡羅(AKMC)、實體動力學(xué)蒙特卡羅(OKMC)和事件動力學(xué)蒙特卡羅(EKMC)。在微結(jié)構(gòu)演化模擬中，OKMC方法對原子運動進行粗?；?，不考慮原子位置分布，將缺陷視作在晶格中移動的實體，用抽取事件發(fā)生的形式來模擬缺陷的演化，提高模擬尺度。本文將對OKMC算法和模型進行介紹，通過一些典型例子來介紹OKMC在材料輻照模擬方面的應(yīng)用，從這些例子中分析OKMC方法的優(yōu)越性和有待提高發(fā)展的方面。

1 OKMC的算法和模型介紹

通常，OKMC方法將缺陷及缺陷簇看作一個具有捕獲半徑的實體，實體具有與之關(guān)聯(lián)的事件及事件可能發(fā)生的概率，通過事件的概率來抽取事件的發(fā)生，推進相應(yīng)時間步長來模擬微結(jié)構(gòu)演化。輻照導(dǎo)致微結(jié)構(gòu)演化的多尺度特性示于圖1。

圖1 輻照導(dǎo)致微結(jié)構(gòu)演化的多尺度特性Fig.1 Multiscale characteristics of microstructure evolution under irradiation

1.1 KMC算法簡介

不同于經(jīng)典的基于Metropolis算法[4]的MC方法，在輻照領(lǐng)域應(yīng)用的KMC方法多是無拒絕的MC方法，多基于RTA(剩余時間算法，也稱為BKL算法)[5]。對于一個給定的系統(tǒng)，根據(jù)事件發(fā)生的速率所占權(quán)重來執(zhí)行事件抽取，完成狀態(tài)轉(zhuǎn)化并計算間隔時間。在KMC中事件i發(fā)生的速率Γi可由下式計算：

(1)

式中：ΔEi為事件發(fā)生的勢壘；kB為玻爾茲曼常數(shù)；T為溫度；vi為嘗試頻率，可根據(jù)諧波過渡態(tài)理論得出：

(2)

時間步長δt和平均時間步長Δt可由下式計算：

δt=-lnr/Γtot

(3)

Δt=1/Γtot

(4)

式中：Γtot為所有事件反應(yīng)速率之和；r為0到1的隨機數(shù)。

KMC將“原子運動軌跡”粗?；癁椤绑w系組態(tài)躍遷”，忽略了若干次不會引發(fā)組態(tài)躍遷的原子振動，給出單位時間內(nèi)體系發(fā)生躍遷的概率，從而推導(dǎo)出發(fā)生1次組態(tài)躍遷的時間。因此KMC的時間步長可以看作真實物理系統(tǒng)中體系穿越勢壘發(fā)生組態(tài)變化所需要的若干次振動的總時間。

因此，KMC算法的基本步驟可總結(jié)如圖2所示。

圖2 KMC算法流程圖Fig.2 Flowchart of KMC algorithm

1.2 OKMC模型

OKMC模型主要包括實體以及實體之間可能發(fā)生的反應(yīng)，諸如自間隙原子、雜質(zhì)原子、嬗變元素、空位以及它們形成的團簇、位錯、空洞等缺陷均被當成實體。OKMC實時跟蹤這些實體(缺陷)的位置。實體的特征包括它們的類型、尺寸、位置、方向、形狀和捕獲半徑。OKMC方法可將實體處理成球形或其他任意形狀[6-8]。另外，在OKMC方法中，晶界和表面被視作面吸收缺陷，能夠吸收到達其表面的其他缺陷[6,9-10]。實體的捕獲半徑通常使用第一性原理進行計算[11]，當實體之間的距離小于實體的捕獲半徑之和時，實現(xiàn)實體之間的反應(yīng)，如成團、復(fù)合或湮滅。

OKMC中的事件由用戶自定義。通常，事件庫中包含的事件如下。

1) 擴散(遷移)：即缺陷向近鄰位置的移動。根據(jù)材料的不同、缺陷的不同，擴散會有不同的方向、不同的維度，這些擴散事件的發(fā)生由相應(yīng)的速率控制。

2) 轉(zhuǎn)向：一維或二維運動的缺陷改變當前擴散方向。

3) 成團：相同種缺陷或不同種缺陷之間的距離小于它們的捕獲半徑之和時，形成一個更大的缺陷團簇，如兩個單空位(V)發(fā)生成團形成1個含有兩個空位的空位團簇，即V1+V1→V2。

4) 復(fù)合：異種缺陷之間的合并。如1個三空位的團簇和1個單自間隙原子(SIA)發(fā)生復(fù)合之后，形成1個雙空位團簇，即V3+SIA1→V2。

5) 解離：缺陷團簇釋放出單個缺陷的過程。如1個三空位團簇釋放出1個單空位的過程，V3→V2+V1。

6) 捕獲：缺陷被一些特殊的雜質(zhì)原子或位錯等捕獲。

7) 湮滅：指缺陷被自由表面、晶界或其他微結(jié)構(gòu)吸收移除的過程。

8) 陷阱突變：金屬中含有一些稀有氣體元素時觀察到的常見現(xiàn)象。如在聚變和裂變條件下，金屬材料中含有氦時，氦原子被金屬原子排斥，形成團簇。當團簇中氦原子的數(shù)量較多時，氦團簇會擠出1個或多個基質(zhì)原子，產(chǎn)生1個或多個弗倫克爾對(Frenkel Pair)，甚至產(chǎn)生位錯環(huán)。

9) 輻照事件：在OKMC中通過不斷引入缺陷模擬輻照過程，可通過隨機引入缺陷的方式實現(xiàn)，也可通過以一定速率引入缺陷級聯(lián)的方式實現(xiàn)。引入速率可根據(jù)NRT模型[12]計算，即不同反沖能量T下的單個級聯(lián)產(chǎn)生的離位次數(shù)vNRT為：

(5)

式中，ED為材料原子離位閾能。

其中缺陷的成團和復(fù)合屬于自發(fā)事件，當缺陷之間的捕獲半徑之和小于缺陷之間的距離時，便會發(fā)生。而缺陷的擴散、轉(zhuǎn)向、解離是由一系列的勢壘和嘗試頻率來描述的，其中擴散事件發(fā)生的速率為：

(6)

式中：vd為擴散事件的嘗試頻率，金屬中，單個點缺陷遷移的嘗試頻率約為金屬的德拜頻率；Ed為相應(yīng)的擴散勢壘，通常由第一性原理[11,13]和分子動力學(xué)計算[14]可得出。

轉(zhuǎn)向事件發(fā)生的速率為：

(7)

式中：vt為轉(zhuǎn)向事件的嘗試頻率；Et為相應(yīng)的轉(zhuǎn)向勢壘。

解離事件發(fā)生的速率為：

(8)

式中：ve為解離事件的嘗試頻率；Eb為相應(yīng)缺陷的結(jié)合能。

2 OKMC方法在材料輻照領(lǐng)域中應(yīng)用的典型例子

OKMC方法在輻照模擬領(lǐng)域應(yīng)用甚廣，既可為更高尺度的模擬方法提供輸入，也可在大尺度范圍(空間尺度可達微米級別，時間尺度可達月、年)模擬材料的輻照效應(yīng)。本節(jié)按照OKMC方法在不同尺度下的應(yīng)用進行介紹。

2.1 計算缺陷阱吸收強度

OKMC方法可計算缺陷阱吸收強度，來為平均場速率理論這類更高尺度的模擬方法提供必要參數(shù)，因此獲得可靠的缺陷阱吸收強度對于平均場速率理論的計算具有重大的意義。在速率理論中，缺陷阱吸收強度指可移動的缺陷被一些具有特定形狀、大小、可充當陷阱的缺陷團簇、微結(jié)構(gòu)吸收的速率[15-18]。缺陷阱吸收強度與缺陷的平均自由程的平方呈反比。它與缺陷的大小、形狀和密度以及缺陷遷移的維度有關(guān)?；跀U散理論，Brailsford和Bullough推導(dǎo)了三維極限遷移下的球形缺陷[15]和晶界[16]的缺陷阱吸收強度理論值。Barashev等[19]推導(dǎo)了在一維極限遷移下缺陷的缺陷阱吸收強度理論值。Wiedersich[20]發(fā)表了關(guān)于位錯線的缺陷阱吸收強度的研究結(jié)果，隨后Nichols[21]對此進行了修正。在一維到三維的過渡遷移機制中，缺陷阱吸收強度的理論推導(dǎo)值可由Trinkaus和Heinisch等[22-25]提出的主曲線函數(shù)進行描述。

一般來說，在使用OKMC方法模擬計算缺陷阱吸收強度時，模擬盒子中1次僅放入1個可移動的缺陷，直到這個可移動缺陷被事先預(yù)置的吸收體捕獲之后，再引入下一個缺陷并記錄該缺陷跳躍的次數(shù)。經(jīng)過足夠多次的模擬之后，得到的缺陷阱吸收強度k2如下式所示：

(9)

式中：n為缺陷擴散的維度；d為缺陷跳躍1次的距離；〈N〉為缺陷的平均跳躍次數(shù)。

根據(jù)經(jīng)驗，要使獲得的缺陷阱吸收強度的誤差小于1%，至少需統(tǒng)計104次，要使誤差小于0.1%，則需要統(tǒng)計106次以上[26]。在OKMC中，一維遷移的缺陷始終在一個方向上來回進行跳躍，三維遷移的缺陷在每跳躍一步之后，可以隨機選擇跳躍方向。OKMC處理一維到三維的過渡遷移機制主要有兩種方法：一是在遷移固定的次數(shù)nd之后使缺陷改變跳躍方向；二是給缺陷分配轉(zhuǎn)向能壘Er，讓缺陷在每跳躍一步之后有一定的概率改變跳躍方向。

Jansson等[9]使用OKMC方法計算了位錯和位錯線的缺陷阱吸收強度，并且與理論值符合良好。Malerba等[27]則使用OKMC方法分別計算了球形吸收體、晶界對不同遷移維度的缺陷的缺陷阱吸收強度。一維和三維運動的情況下可直接與理論公式進行對比，混合一維/三維運動情況下，需與主曲線方程進行對比。對于球形吸收體和球形晶界，OKMC模擬計算得出的缺陷阱吸收強度均與理論值符合較好。

另外，有研究[22-23,27-28]表明：OKMC的仿真結(jié)果會高估大體積吸收體的缺陷阱吸收強度，而低估小體積吸收體的缺陷阱吸收強度；與大體積分數(shù)情況相比，在小體積分數(shù)下，周期性分布的位錯線缺陷阱吸收強度與理論值吻合更好。針對這些偏差，Hou等[29]做了相關(guān)的研究，揭示了其間的差異起源，并且在理論上提出了修正。他們認為，低體積分數(shù)的情況下，造成偏差的主要原因是：OKMC模型中離散的擴散形式與實際連續(xù)擴散之間的差異；高體積分數(shù)的情況下，吸收體之間的相互重疊情況是造成差異的主要原因；另外，OKMC模型中位錯線的周期性分布是導(dǎo)致一維遷移下出現(xiàn)鋸齒狀差異的原因。在考慮了這些差異并做出修正之后，模擬值與理論值之間得到了更好的吻合。Ahlgren等[26]對缺陷阱吸收強度的理論值推導(dǎo)進行了修正，并開發(fā)了名為“N-jumps”的蒙特卡羅方法：通過可移動缺陷與吸收體之間的最小距離來確定可能需要的跳躍次數(shù)，然后使用1次“大的跳躍”等效替代N次單獨跳躍的效果，從而提升模擬速度。在吸收體的體積分數(shù)較低時，提升效果更為顯著。

通過以上的研究可了解到，OKMC方法是模擬計算缺陷阱吸收強度的有效工具，并且與理論推導(dǎo)值有很好的吻合。但在某些吸收體體積分數(shù)下，OKMC方法的仿真結(jié)果與理論值總是存在一些偏差，這種偏差來源于OKMC模型的離散化處理與理論連續(xù)性假設(shè)之間不可避免的差異。因此，有效的處理方法是分析差異來源，并且基于現(xiàn)有的模型進行適當修正。OKMC方法缺陷阱吸收強度計算的精準度取決于對缺陷動力學(xué)模型細致的考量，因此需要MD或DFT提供更為準確的缺陷互相作用模型，以及缺陷的熱力學(xué)和動力學(xué)參數(shù)。此外，在小吸收體體積分數(shù)和缺陷一維遷移條件下，缺陷被吸收前的平均自由程較長，導(dǎo)致OKMC模擬緩慢。如何有效地解決該問題，也是OKMC方法需要改善的方向之一。

2.2 級聯(lián)內(nèi)退火模擬

輻照產(chǎn)生的PKA與其他原子相互碰撞，導(dǎo)致材料內(nèi)部原子發(fā)生離位和進一步的碰撞，形成一定的缺陷空間分布，這個過程稱為位移級聯(lián)。隨著級聯(lián)能量的耗散，產(chǎn)生的缺陷集中分布在初級損傷區(qū)域，這個過程通常在十幾皮秒內(nèi)完成。隨后，初級損傷區(qū)域內(nèi)的缺陷之間會發(fā)生進一步的復(fù)合、成團等相互作用，直到缺陷的空間、尺寸分布趨于穩(wěn)定，這個過程稱為級聯(lián)內(nèi)退火。級聯(lián)內(nèi)退火過程的持續(xù)時間通常為數(shù)納秒到數(shù)十納秒不等。位移級聯(lián)經(jīng)過級聯(lián)內(nèi)退火之后尚存的缺陷數(shù)目、大小和空間分布以及材料中的雜質(zhì)會對輻照微結(jié)構(gòu)演化產(chǎn)生影響[30]。因此，充分了解級聯(lián)內(nèi)退火的特征對級聯(lián)插入式的OKMC輻照模擬具有重大意義。

Xu等[31]使用分子動力學(xué)模擬得到的初級輻照損傷數(shù)據(jù)庫進行了bcc鐵中的級聯(lián)退火模擬，并進行了參數(shù)敏感度研究。得到的結(jié)果可作為平均場速率理論的源項。模擬主要結(jié)果如圖3所示，其中自由的SIA分數(shù)表示經(jīng)過級聯(lián)之間的相互作用之后可脫離級聯(lián)區(qū)域并且不再與級聯(lián)內(nèi)其他缺陷發(fā)生反應(yīng)的SIA分數(shù)?？煽闯鯬KA能量越低，模擬結(jié)果的波動越大。這是由于低PKA能量下，級聯(lián)中的缺陷數(shù)量較少，在統(tǒng)計上的波動更大。而在能量大于20 keV時，統(tǒng)計誤差就較小了。從結(jié)果可看出，自由的SIA分數(shù)約為55%～70%，與PKA能量的關(guān)聯(lián)度不高。參數(shù)敏感度分析涵蓋了不同的PKA能量、不同的SIA運動機制、缺陷團簇解離情況、不同的捕獲半徑以及模擬盒子的大小。研究發(fā)現(xiàn)，SIA的運動機制對模擬結(jié)果的影響最大，當越多的SIA團簇以三維形式擴散時，逃離級聯(lián)的SIA比例越低。另外，當所有的SIA團簇均以三維形式擴散時，復(fù)合效率最高。

圖3 bcc鐵中不同PKA能量的級聯(lián)退火模擬結(jié)果[31]Fig.3 Cascade annealing simulation results of different PKA energy in bcc iron[31]

Nandipati等[6,32]使用OKMC對鎢的級聯(lián)內(nèi)退火進行了模擬，并且進行了點缺陷動力學(xué)參數(shù)敏感度的分析。同樣以分子動力學(xué)模擬[33]結(jié)果作為輸入，典型模擬結(jié)果如圖4所示。級聯(lián)退火過程可分為兩個階段：階段1，SIA和空位之間的復(fù)合基本完成；階段2，SIA逐漸擴散到吸收邊界而從模擬盒子中移除，而空位基本保持不動，因此模擬盒子中SIA數(shù)量不斷減少，而空位數(shù)量基本不變。同時，Nandipati等也進行了參數(shù)敏感度研究，發(fā)現(xiàn)級聯(lián)內(nèi)退火行為取決于缺陷的空間和尺寸分布，以及缺陷的相對遷移速率。具體來說：1) PKA能量會影響缺陷的尺寸分布和空間分布，進而影響復(fù)合效率。2) 復(fù)合效率隨溫度從300 K到1 025 K呈現(xiàn)增長的趨勢，然后在2 050 K時下降，與溫度的關(guān)系表現(xiàn)出倒U型數(shù)曲線，這說明退火行為的演化是溫度和PKA能量共同作用的結(jié)果。3) 在金屬鎢中，由于SIA的遷移速率遠大于空位，因此SIA的遷移和小SIA團簇的轉(zhuǎn)向主導(dǎo)著級聯(lián)內(nèi)退火行為；SIA擴散的維度也存在影響，三維擴散下的退火效率高于一維擴散。

圖4 鎢中PKA能量為75 keV的級聯(lián)在300 K的溫度下退火10 ns點缺陷的數(shù)量變化[6]Fig.4 Change of number of point defects for cascades of PKA energy of 75 keV in tungsten annealed at temperature 300 K for 10 ns[6]

級聯(lián)內(nèi)退火模擬服務(wù)于相應(yīng)條件下材料微結(jié)構(gòu)演化模擬，模擬結(jié)果經(jīng)處理也可作為平均場速率理論的關(guān)鍵輸入。但由于單級聯(lián)退火的時間通常很短，一些在單級聯(lián)退火過程中不敏感的參數(shù)，在長期演化中仍可能會對模擬造成不可忽略的影響。因此，如何有效地將單級聯(lián)退火的模擬結(jié)果應(yīng)用到長時間演化模擬中，需進行進一步的研究。另外由于OKMC方法忽略了缺陷的原子尺度的細節(jié)，缺少了缺陷構(gòu)型的信息，而不同構(gòu)型缺陷的動力學(xué)參數(shù)往往也不相同。恰是由于OKMC方法在這方面的局限性，所以在模擬中不得不進行一些假設(shè)。

2.3 微結(jié)構(gòu)演化模擬

級聯(lián)碰撞產(chǎn)生的缺陷的空間分布和遷移率是不均勻的，這種不均勻性導(dǎo)致局部的缺陷之間發(fā)生相互作用，這種行為被稱為空間關(guān)聯(lián)效應(yīng)[34]?？臻g關(guān)聯(lián)效應(yīng)的存在，會使缺陷的空間、尺寸分布變復(fù)雜，對長時間的微結(jié)構(gòu)演化模擬產(chǎn)生很大的影響。不同于一般的平均場方法，OKMC方法可在考慮到空間關(guān)聯(lián)性的前提下，模擬空間尺度為納米到微米級別、時間尺度根據(jù)不用的應(yīng)用場景可到小時、日，甚至月、年的微結(jié)構(gòu)演化。模擬可得到材料內(nèi)部缺陷的尺寸分布和空間分布，并且能夠與實驗結(jié)果進行對比驗證。本節(jié)將通過反應(yīng)堆壓力容器鋼的脆化機理研究、碳雜質(zhì)對微結(jié)構(gòu)演化的影響、空洞陣列的形成3例，對OKMC方法在模擬微結(jié)構(gòu)演化上的應(yīng)用進行介紹。

1) 反應(yīng)堆壓力容器鋼的脆化機理研究

用于輻照實驗和反應(yīng)堆運行的材料很難做到100%的純度，在輻照過程中也會有嬗變元素的生成，加之合金材料的廣泛應(yīng)用，在OKMC方法中如何處理多種元素成為了研究的熱點，如反應(yīng)堆壓力容器鋼的脆化機理研究[35-40]、氫元素和氦元素對鎢微結(jié)構(gòu)演化的影響[41-42]、碳雜質(zhì)對空洞腫脹的影響[36,43-44]等。

目前對于異間隙原子(FIA)的處理主要分為顯式引入和隱式引入兩種方式。前者通過直接在模型引入原子，增加相應(yīng)的FIA與點缺陷形成的絡(luò)合物和相互作用；后者通過等效替代或引入缺陷阱吸收強度的方法來實現(xiàn)。

反應(yīng)堆壓力容器(RPV)作為核電站反應(yīng)堆的重要組件，其安全性直接影響著反應(yīng)堆的使用壽命。通常使用低合金貝氏體鋼、鐵素體鋼作為RPV的結(jié)構(gòu)材料。長時間暴露于強輻射下，RPV鋼中會產(chǎn)生大量的點缺陷和小缺陷簇，而這些缺陷簇可能演變成空洞和位錯環(huán)，進而使RPV在宏觀上表現(xiàn)為屈服強度增加(硬化)、斷裂韌性降低、韌脆轉(zhuǎn)變溫度(DBTT)增加(脆化)[45-50]。

RPV鋼的脆化主要由于富錳(Mn)和富鎳(Ni)團簇的形成[51-52]。OKMC方法由于具有大時間尺度的特性，足以模擬RPV鋼的整個壽命周期，Chiapetto等[35,37-38,40]對此提出了Fe-C-MnNi合金模型，使用OKMC方法研究了溶質(zhì)元素在輻照下對納米級微結(jié)構(gòu)演化的影響，通過修改缺陷遷移參數(shù)隱式地引入合金元素。模擬結(jié)果表明，與Fe-C相比，F(xiàn)e-C-MnNi合金中大的空位團簇、SIA團簇和空洞的形成受到抑制。這些結(jié)果與正電子湮沒譜(PAS)實驗值一致[53]。此外，Mn濃度也會對微結(jié)構(gòu)演化產(chǎn)生影響。Mn的濃度升高會導(dǎo)致SIA團簇增多，但在Mn的濃度高于1.2%之后，SIA團簇的數(shù)量不再明顯受到Mn濃度變化的影響。另外，由于空位團簇的遷移率不隨Mn濃度改變而變化，所以并未觀察到空位團簇的明顯變化。

Fe-C-MnNi的OKMC模擬可清楚地揭示RPV鋼脆化的機制：點缺陷簇遷移率的降低，導(dǎo)致Fe-C和Fe-C-MnNi合金之間微結(jié)構(gòu)演化發(fā)生顯著差異。為更好地與實驗結(jié)果進行比對，Chiapetto等[35]對這套模型和部分參數(shù)集進行了進一步的修正，然后針對兩種不同的合金(Fe-C-MnNi合金和Ringhals RPV鋼)進行了模擬，模擬結(jié)果與實驗結(jié)果符合良好，如圖5所示。此研究進一步證實了富MnNi團簇是由小間隙原子環(huán)形成的溶質(zhì)富集產(chǎn)生的。

圖5 Fe-C-MnNi合金和Ringhals RPV鋼中SIA數(shù)密度隨輻照劑量的變化趨勢[35]Fig.5 Density of SIA in Fe-C-MnNi alloys and Ringhals RPV steels with irradiation dose[35]

2) C雜質(zhì)對微結(jié)構(gòu)演化的影響

C雜質(zhì)廣泛存在于各種金屬中，無論是在RPV鋼中，還是在其他金屬材料，如選用為ITER偏濾器材料的金屬鎢[54]中，C雜質(zhì)的存在均會對其機械性能造成影響[45,51]。

OKMC模型的隱式引入方法中，使用球形吸收體來表示等效替代C原子效應(yīng)的陷阱元素。在α-Fe中，C原子易與空位結(jié)合形成CV、CV1、CV2、C2V等絡(luò)合物，影響SIA和空位在基質(zhì)中的遷移。相關(guān)的MD模擬發(fā)現(xiàn)CV、CV2和C2V絡(luò)合物非常穩(wěn)定，并且，CV絡(luò)合物與1/2〈111〉SIA簇的結(jié)合能取決于與SIA簇的作用位置是邊緣還是中心[55]。而在不同的退火溫度下，絡(luò)合物的形式也有不同[36]。在α-Fe中，Jansson等[36,44]使用隱式引入的方法來模擬C原子的影響。首先通過簡單的預(yù)實驗，顯式引入C原子，確定CV絡(luò)合物在不同溫度下的主要構(gòu)型，然后在不同的溫度下確定不同尺寸的SIA團簇與C、CV絡(luò)合物的結(jié)合能。使用陷阱元素等效替代C以及CV絡(luò)合物之后，陷阱與點缺陷之間的捕獲效應(yīng)體現(xiàn)在缺陷團簇的動力學(xué)參數(shù)變化上。點缺陷一旦被陷阱捕獲，便以結(jié)合能Eb來束縛。

圖6 空位團簇數(shù)密度和平均尺寸與輻照劑量的關(guān)系[44] Fig.6 Relationship between density and average size of vacancy clusters and irradiation dose[44]

圖6為Jansson等模擬α-Fe中Fe-C的微結(jié)構(gòu)演化的部分結(jié)果。模擬溫度為343 K，采用級聯(lián)插入的方式模擬輻照，缺陷級聯(lián)來自于已有數(shù)據(jù)[56-61]，引入速率相當于7×10-7dpa/s?？瘴粓F簇的數(shù)密度隨輻照劑量的上升而上升，與實驗結(jié)果符合良好。同時，Jansson等還進行了等時退火模擬，起始溫度為333 K，溫度梯度為50 K，在一個溫度下退火3 600 s。依據(jù)預(yù)實驗的結(jié)果，改變不同溫度下C源陷阱與SIA的結(jié)合能?？傮w上，模擬結(jié)果與實驗值符合較好。

OKMC的靈活性可滿足不同方式的隱式引入方法。Castin等[10]通過計算在不同濃度下C的缺陷阱吸收強度，然后將其轉(zhuǎn)化為點缺陷與C源陷阱的作用概率，來表示C雜質(zhì)引入的效果。并且，此研究還假設(shè)：隨空位團簇的長大，C對空位團簇的捕獲效應(yīng)降低，即相應(yīng)的結(jié)合能Eb降低。研究結(jié)果表明，C含量的變化對空洞形成的動力學(xué)過程有很大影響。因為C對空位的捕獲導(dǎo)致了進入空洞的空位減少，使空洞平均直徑減小，延緩了空洞的生長。另外，由于空位成團的減少，小空位簇的彌散程度升高，導(dǎo)致輻照引入的SIA易與之復(fù)合，使空洞數(shù)密度降低。

以上OKMC研究表明了雜質(zhì)元素通過與點缺陷的作用來影響點缺陷團簇的演化，從而影響微結(jié)構(gòu)的演化。但這些OKMC研究大都對模型進行了一些假設(shè)或簡化處理，如假設(shè)SIA團簇與C的作用存在尺寸閾值，這個閾值決定了二者的作用類型是邊緣弱作用還是中心強作用[44]。這種假設(shè)雖是根據(jù)理論研究結(jié)果和DFT模擬結(jié)果進行的簡化處理，但其正確性尚有待驗證，也需進行相關(guān)的敏感度研究予以支撐。

在使用OKMC方法處理溶質(zhì)元素的影響時，顯式地引入新的元素能夠直接處理溶質(zhì)元素在基體中的擴散，得到直觀的、明確的團簇尺寸和空間分布。但實際上很難實現(xiàn)溶質(zhì)原子與空位或SIA缺陷之間相互作用的精確建模，并且這樣的處理給計算增加了負擔。首先，根據(jù)KMC算法可知，顯式地將它們包括在模型中，會大幅增加總反應(yīng)速率Γtot，減小時間步長Δt，將不可避免地并顯著增長KMC總模擬時間。其次，模型中的事件數(shù)量將呈現(xiàn)幾何式的增長，變得難以管理。隱式引入方法使得模型變得更易處理，減輕了計算負擔，但打破了溶質(zhì)元素在基體中分布的不均勻性，不能利用OKMC處理空間關(guān)聯(lián)性的優(yōu)勢，對模擬結(jié)果的影響程度尚有待研究。因此，在保留空間關(guān)聯(lián)性的前提下，如何融合其他方法來提高OKMC模擬效率也是有價值的研究方向之一。

3) 空洞陣列模擬

空洞陣列是指受輻照材料在高溫、高劑量輻照下形成的空洞有序排列的情況，空洞超晶格的構(gòu)型通常與材料本身的晶格構(gòu)型相同。1971年，Evans在鉬離子輻照實驗中[62]首次發(fā)現(xiàn)空洞陣列的形成。1972年，Sikka等發(fā)表了鎢在中子輻照下形成空洞點陣的實驗研究結(jié)果[63]。這項研究在EBR-Ⅱ反應(yīng)堆中進行，于550 ℃照射至1×1022cm-2中子通量。他們分析得到了空洞超晶格參數(shù)，并且與金屬Mo和Ta進行了比較。Loomis等在金屬Nb及其合金中研究了空洞點陣的晶格參數(shù)隨輻照條件的變化[64]。Tanno等[65]進行的鎢中子輻照實驗中，在750 ℃輻照至1.54 dpa時，也發(fā)現(xiàn)了空洞陣列的形成。

目前關(guān)于空洞陣列形成的主流觀點是陰影效應(yīng)[24]及其衍生理論。陰影效應(yīng)認為，形成空洞陣列的核心點在于空位和SIA運動方式的不同。如在bcc金屬鎢中，空位的擴散一般認為是三維的，而可移動的SIA沿著密排方向一維遷移。恰是由于這種運動方式的差異使得空位在特定位置聚集、排列，最終形成空洞陣列。

OKMC模擬能夠很好地體現(xiàn)SIA不同運動方式的差異，有效地處理空洞形核過程，并且能夠給出空洞的尺寸、空間分布，很適合用于模擬空洞陣列的形成。圖7示出了使用OKMC方法模擬空洞陣列的結(jié)果示意圖。Heinisch等[24]使用OKMC方法針對空洞陣列的形成機理進行了研究。他們設(shè)計了1組對照模擬實驗：預(yù)置空洞并按陣列排布、預(yù)置空洞隨機分布。模擬盒子中僅入射SIA，發(fā)現(xiàn)空洞陣列的形成需要足夠長的SIA一維運動路徑，而SIA過于頻繁的轉(zhuǎn)向?qū)⒉粫纬煽斩搓嚵小Ａ?組對照模擬實驗分別在3種方案下進行模擬：1) 預(yù)置14個按fcc結(jié)構(gòu)排布的大空洞作為“種子”；2) 在模擬盒子中心預(yù)置1個大空洞作為“種子”；3) 無預(yù)置空洞作為“種子”。3種方案均預(yù)先隨機放置了1 000個小空洞，在模擬盒子中持續(xù)入射空位和SIA來模擬輻照。以上3種方案均形成了空洞陣列，但方案1中形成的空洞陣列的晶格參數(shù)最大，并且在有“種子”的情況下，空洞會發(fā)生解離或移位到格點上。該研究側(cè)重于探索在陰影效應(yīng)的理論基礎(chǔ)下，SIA的運動方式如何對空洞陣列形成穩(wěn)定性造成影響。模擬中采用了許多假設(shè)，如入射的SIA和空位均視作尺寸為7的團簇、缺陷遍歷模擬盒子4次便會被刪除等，并且忽略了空洞的形核過程。這與現(xiàn)今許多OKMC模擬一致，并未采用真實的模型和參數(shù)，缺少與實驗之間的對比來驗證理論的正確性。

圖7 OKMC方法模擬空洞陣列的形成Fig.7 Formation of void lattice simulated by OKMC method

因此，目前關(guān)于空洞陣列的OKMC模擬尚有較大的發(fā)展空間，特別是在模擬盒子形狀和尺寸、邊界條件的設(shè)置問題上。直接以真實的晶粒尺寸作為模擬盒子的尺寸，使用吸收性邊界條件作為晶界是模擬多晶材料輻照最理想的設(shè)置，但這樣會導(dǎo)致計算量劇增，不太現(xiàn)實。所以，周期性邊界條件受到廣泛應(yīng)用。然而，OKMC模型在處理一維遷移的缺陷時，如果使用周期性邊界條件，在缺陷穿越邊界后人為地限定重新進入模擬盒子的位置，將不可避免地導(dǎo)致一維運動缺陷的遷移路徑受到模擬盒子邊長比例的影響。如果使用立方體的模擬盒子，會使一維運動的缺陷一直重復(fù)同一條路徑，如圖8a所示；如果使用非立方體的模擬盒子，如圖8b、c所示，會使一維遷移的缺陷每次穿越邊界之后的路徑之間的間距與模擬盒子邊長的最大公約數(shù)相等。這樣會使得OKMC模擬中形成的空洞陣列的晶格參數(shù)將會被模擬盒子邊長比所主導(dǎo)，甚至?xí)谏w輻照條件和缺陷動力學(xué)參數(shù)的影響。因此，如何合理地去除周期性邊界條件對模擬一維運動的影響也是OKMC方法在材料輻照模擬領(lǐng)域的發(fā)展方向之一。

圖8 不同邊長比的情況下一維運動缺陷的軌跡Fig.8 Trajectory of one-dimensional motion defects in simulation boxes with different length-to-width ratios

3 總結(jié)與展望

本文對OKMC方法的原理和建模進行了介紹，通過例舉一些現(xiàn)有研究，對OKMC方法在材料輻照領(lǐng)域的應(yīng)用進行了詳細介紹，并對OKMC方法的優(yōu)勢以及待發(fā)展和解決的問題進行了闡釋?？偟膩碚f，OKMC方法具有MD、CD等方法所無法替代的特性和優(yōu)勢，主要體現(xiàn)在以下幾點。

1) OKMC模型除具有足夠的時空尺度之外，還具備捕獲缺陷空間關(guān)聯(lián)性的能力，而這些空間關(guān)聯(lián)性極有可能是影響受輻照材料微結(jié)構(gòu)演化的關(guān)鍵因素。

2) OKMC模擬能夠給出可與實驗進行直接對比的缺陷尺寸、空間分布等信息，并且現(xiàn)有的研究表明，OKMC模擬結(jié)果能夠與實驗值符合較好。

3) OKMC可用于確定平均場理論模型中的源項或輸入項。

4) OKMC模型靈活，可建立缺陷參數(shù)與輻照效應(yīng)之間的橋梁，很適合用于研究參數(shù)、假設(shè)對損傷退火和損傷累積的影響。

但OKMC仍存在不足。首先，由于OKMC方法不具備預(yù)知性，其誤差主要來源于事件庫和速率的不完整性，其次是低尺度輸入的不確定性以及模型的統(tǒng)計誤差。在真實的系統(tǒng)中，事件發(fā)生的速率是會隨著時間變化的，并且還與系統(tǒng)的狀態(tài)相關(guān)，而在OKMC中該速率與時間和狀態(tài)無關(guān)。另外，OKMC中還假設(shè)各事件之間是獨立的，然而在真實系統(tǒng)中，它們可能是相關(guān)的，這也是誤差的來源之一。但完備的事件庫往往不易獲得，尤其是當所研究的體系性質(zhì)復(fù)雜并且涉及多種組元或復(fù)雜的晶體結(jié)構(gòu)時。

因此，為滿足模擬需求，需對模型進行必要的簡化，在準確性與效率中尋找平衡。在這種情況下對缺陷動力學(xué)模型進行合理假設(shè)，并且開展相應(yīng)的敏感性分析是很有必要的。如Nandipati等[32]展開的敏感性分析表明了級聯(lián)的退火效率主要受SIA遷移的影響，其中小SIA團簇的轉(zhuǎn)向和遷移顯著影響著SIA與空位的復(fù)合。因此，在建模過程中對于SIA的遷移和轉(zhuǎn)向需要詳細刻畫，增加模擬的準確性；而對于其他的事件，可做一些簡化，提高模擬的效率。此外，也可通過與一些簡單實驗進行對比，反推出模型中刻畫不準確的參數(shù)和事件，調(diào)整相應(yīng)的參數(shù)值和事件庫。通過反復(fù)迭代進行對比調(diào)整，校正事件庫和參數(shù)，增強OKMC方法在復(fù)雜條件下的模擬準確性。

即時動力學(xué)蒙特卡羅(On-the-fly KMC)算法[66-69]可克服事件庫完備性的問題，即在執(zhí)行KMC計算之前無需給出可能的事件及其發(fā)生速率，而是在算法步驟執(zhí)行的同時確定可能的事件并計算其速率。這種方法在模擬合金中缺陷的演化尤其有效。因為合金中的一些特定的反應(yīng)速率，特別是缺陷的擴散速率，取決于局部環(huán)境。但這種算法在OKMC方法中難以實現(xiàn)，通常是在AKMC方法中進行應(yīng)用。因為實時確定可能的事件并計算其速率需知體系中所有原子的位置，以及它們之間的互相作用勢這些信息，但OKMC方法中所有缺陷均被視為實體，丟失了這些原子尺度的信息。因此，未來的研究可考慮將OKMC和具備即時算法的AKMC結(jié)合，通過AKMC方法來實時更新關(guān)鍵的參數(shù)值和可能發(fā)生的事件，從而達到完備OKMC的事件庫、提高OKMC模擬準確性的目的。

另外，由于OKMC中時間步長與總速率呈反比，因此在高溫和高劑量率的條件下，OKMC方法的計算步驟量將會異常龐大，所能模擬的總輻照劑量將很小，不能滿足與實驗對照的需求。因此，這種情況下進行OKMC算法的并行化就十分有必要了[70-71]。但在處理一些高度不均勻的系統(tǒng)或系統(tǒng)內(nèi)不同缺陷的遷移速率相差甚大的情況時，并行算法的性能將會嚴重受限。

綜上所述，未來OKMC方法在輻照材料模擬領(lǐng)域應(yīng)用需克服的挑戰(zhàn)主要包括：

1) 解決算法的并行化的難題，顯著提高OKMC方法的模擬效率；

2) 完善OKMC模型以處理多組分材料系統(tǒng)，包括不斷發(fā)生嬗變的系統(tǒng)；

3) 與不同時間尺度、空間尺度的方法進行結(jié)合，構(gòu)建材料輻照模擬的完整平臺。