尚子豪，商紅慧，王東杰，張云泉，賀新福，陳澤華，王 棟，張廣婷

（1.太原理工大學(xué)大數(shù)據(jù)學(xué)院，山西 太原030024；

2.中國科學(xué)院計算技術(shù)研究所計算機體系結(jié)構(gòu)國家重點實驗室，北京100190；

3.中國原子能科學(xué)研究院，北京102413）

1 引言

核能是一種可以大規(guī)模使用的工業(yè)能源，核能發(fā)電具有發(fā)電成本較低、發(fā)電小時數(shù)長、電網(wǎng)接入穩(wěn)定以及零碳排放等優(yōu)勢。核反應(yīng)堆壓力容器RPV（Reactor Pressure Vessel）是防止放射性物質(zhì)泄漏的最主要屏障，并且是反應(yīng)堆中不可拆換的部件，RPV因此成為核電站最關(guān)鍵的部件之一。RPV材料的使用壽命直接決定了核電站的服役壽命，在RPV鋼中加入少量的銅能顯著提高鋼的耐腐蝕性能，同時增加鋼的強度和疲勞抗力等。但是，長時間暴露在高能中子輻照場中會導(dǎo)致RPV鋼中的銅元素以富銅團簇的形式析出，這些高數(shù)量密度的納米富銅團簇是引起RPV鋼產(chǎn)生脆化的主要原因[1，2]。

正常狀態(tài)下，由于銅的自擴散系數(shù)較低，未受到輻照時銅原子基本不能通過擴散聚集而形成團簇，然而高能中子輻照產(chǎn)生的點缺陷會大大提高溶質(zhì)原子的擴散聚集能力，從而產(chǎn)生原子團簇。利用現(xiàn)代實驗技術(shù)，如自擴散示蹤法、電阻率法和正電子湮沒法等[3-5]對點缺陷相關(guān)性質(zhì)的實驗測量，可以推斷出許多有關(guān)缺陷的信息。然而，由于分析方法和儀器分辨率的限制，用這些技術(shù)直接觀察和研究點缺陷團簇的原子結(jié)構(gòu)和性質(zhì)很困難。利用計算機模擬的方法可以對各種點缺陷進行深入的研究。

目前對RPV鋼的富銅團簇的研究主要集中在輻照注量、輻照溫度以及元素含量對其性能的影響，而通過引入點缺陷的方式來觀察富銅團簇析出過程的研究較少。動力學(xué)蒙特卡洛KMC（Kinetic Monte Carlo）方法[6]可以有效地模擬富銅團簇析出的微觀結(jié)構(gòu)演化。許多研究人員致力于開發(fā)不同的KMC方法。然而，這類軟件多為串行版本，這使得對富銅團簇析出過程的模擬只能在較小的體系中進行，也因此只能引入單個或是少量的點缺陷，限制了模擬的精確程度。為此我們自主開發(fā)了可以并行模擬上億原子體系的Open KMC軟件[7]。本文基于Pair勢[8]和嵌入原子勢EAM（Embedded-Atom Method）[9]，通過引入空位點缺陷的方式，模擬熱時效下Fe-Cu二元合金中團簇的析出過程。并對程序中的空位點缺陷選取問題進行了計算優(yōu)化，通過多組測試驗證了優(yōu)化算法的正確性和有效性，同時對優(yōu)化后的程序做了大規(guī)模并行模擬的性能分析。

2 程序優(yōu)化

2.1 基于空位躍遷的KMC方法

KMC模擬方法是一種通過構(gòu)造隨機過程來模擬體系長時間演化的方法。經(jīng)典的原子動力學(xué)蒙特卡洛AKMC（Atomic Kinetic Monte Carlo）方法[10]利用晶格點陣來映射原子與格點的位置關(guān)系，應(yīng)用于提供控制微觀系統(tǒng)演化的精確解。α-Fe的晶格為體心立方（bcc）結(jié)構(gòu)，空位躍遷事件被認為是空位與其周圍第1近鄰原子的交換。圖1分別給出了bcc晶格中心原子的第1、第2和第3近鄰原子（晶格邊長為2）。模擬通常包括4個步驟：（1）計算列表中空位的躍遷概率；（2）計算時間增量；（3）執(zhí)行空位躍遷，增加時間；（4）是否達到時間閾值，若達到則退出，否則循環(huán)執(zhí)行以上步驟。

Figure 1 Three kinds of nearest neighbors for a bcc lattice body center atom圖1 bcc晶格體心原子的三重近鄰

第1步根據(jù)當(dāng)前狀態(tài)計算空位可能發(fā)生的躍遷事件X的概率，這個概率實際上并非真正意義上的0～1的概率，而是一種頻率，其計算方式如下所示：

其中，Γ0代表嘗試頻率，嘗試頻率被設(shè)置為6×1012Hz，k是玻爾茲曼因子，T是絕對溫度，EXa表示事件發(fā)生需要的激活能。考慮到計算時間和結(jié)果精度，我們采用環(huán)境相關(guān)模型[11]，得到了滿足式（2）所示的平衡規(guī)則的激活能EXa：

其中，EX0表示參考激活能，僅取決于與空位交換位置的原子的化學(xué)性質(zhì)。本文中使用的包括Fe：0.65 eV，Cu：0.56 e V。Ei和Ef分別表示空位躍遷前后系統(tǒng)的總能量值。

第2步采用駐留時間算法[12]計算時間增量，它與所有可能發(fā)生的事件頻率之和的倒數(shù)成正比，以便將事件的傳播與每個獨立事件相關(guān)聯(lián)。計算公式如式（3）所示：

其中隨機數(shù)r∈（0，1）用于選擇事件。

第3步是依據(jù)概率隨機確定空位的跳躍方向，發(fā)生跳躍并更新周圍原子的能量和躍遷概率。同時，在模擬時間上累積時間增量Δt。最后一步，重復(fù)上述步驟，直到模擬時間達到預(yù)設(shè)閾值。

在計算第1步中的能量時，原子間相互作用勢采用Pair勢能模型[8]或EAM勢能模型[9]，Pair勢由于其簡潔的形式、較小的計算量和有效的計算結(jié)果，在原子層次的計算機模擬中得到了廣泛的應(yīng)用。Pair勢中原子間的相互作用最遠只考慮到第2近鄰位置。組成系統(tǒng)的所有實體間的化學(xué)相互作用可以通過式（4）描述：

當(dāng)方程中a等于1或2時，分別對應(yīng)于第1或第2近鄰相互作用。Va表示空位，b為Fe-Fe鍵的數(shù)目，c為Va-Va鍵的數(shù)目，l為Fe-Va鍵的數(shù)目，m為Fe-Cu鍵的數(shù)目，n為Va-Cu鍵的數(shù)目，p為Cu-Cu鍵的數(shù)目。

EAM勢中原子間相互作用由Pair勢和多體嵌入勢共同貢獻，并且至多考慮到第3近鄰。其勢能計算方式如式（5）所示：其中，μ和v為原子表示，F(xiàn)（）為多體嵌入能函數(shù)，ρμ為μ處插值的總電子密度。相比于EAM勢能模型，Pair勢能模型的局限性在于缺乏多體嵌入勢能，不能考慮局部環(huán)境的影響。此外，Pair勢能模型需要通過實驗或熱力學(xué)數(shù)據(jù)仔細調(diào)整。本文所用的Pair勢能模型參數(shù)取自文獻[13]，可以得到接近EAM模型的結(jié)果，也可以推廣到更復(fù)雜的多元合金研究當(dāng)中。

2.2 計算優(yōu)化

KMC算法中模型的最主要部分包括躍遷概率和系統(tǒng)能量的計算。如前面提到的，在模擬的第3步空位發(fā)生躍遷之后，需要對該空位周圍原子的能量和躍遷概率進行更新。相當(dāng)于在下次空位躍遷事件發(fā)生之前，我們要對體系內(nèi)所有受到上一次躍遷影響的空位進行更新。以往的KMC版本是將上次發(fā)生躍遷的原子附近理想半徑下的全部空位作為需要被更新的對象。事實上，這些空位當(dāng)中有很多因為距離較遠，而并沒有受到上次躍遷事件的影響，導(dǎo)致模擬過程出現(xiàn)了很多冗余的模擬計算。因此，如何選取需要被更新的空位成為一個關(guān)鍵問題。

被更新空位的選取與采用的勢能函數(shù)的計算方式有關(guān)。對于Pair勢能模型來說，原子間的相互作用最遠考慮到第2近鄰位置。某一空位是否受到上次躍遷影響取決于這個空位的第1近鄰原子是否與上次發(fā)生躍遷的原子存在相互作用。如果能找到發(fā)生躍遷的原子可以影響到的最遠距離的空位位置，并將這個最遠距離作為截斷半徑，每次躍遷事件結(jié)束后只更新截斷半徑以內(nèi)的空位，就能避免大量的冗余計算。

截斷半徑的推導(dǎo)過程在三維坐標(biāo)系中完成。假設(shè)發(fā)生躍遷交換的原子分別為i和j，并在坐標(biāo)系中固定i和j的位置坐標(biāo)。確定Pair勢能模型下最遠處需要被更新的空位位置可以分為下面3個步驟：（1）找到j(luò)原子附近2層近鄰內(nèi)的全部原子并記錄坐標(biāo)；（2）分別計算這些原子與i原子的距離，找到其中一個與i原子距離最遠的原子記為s；（3）記錄s原子的8個1nn近鄰原子坐標(biāo)，并找到其中與i原子距離值最大的對應(yīng)原子記為k。其中近鄰中的nn代表nearest neighbor，表示第幾近鄰的含義。通過以上3步得到最遠處需要被更新的空位k位置。所有原子間位置關(guān)系如圖2所示，將立體圖展開得到圖3中的三視圖，其中晶格邊長為2。

Figure 2 Stereogram of Pair truncation radius derivation圖2 Pair截斷半徑推導(dǎo)立體圖

Figure 3 Three views of Pair truncation radius derivation圖3 Pair截斷半徑推導(dǎo)三視圖

Table 1 Nearest neighbor distance comparison table表1 近鄰距離對照表

通過相同的思路，可以確定EAM勢能模型下空位更新的截斷半徑。與Pair勢能模型的區(qū)別在于，EAM勢能模型的原子間相互作用考慮到第3近鄰位置，因此在尋找空位k步驟的第1步里要找到j(luò)原子附近3層近鄰內(nèi)的全部原子并記錄坐標(biāo)，接著執(zhí)行相同的第2步和第3步。計算出的EAM勢能模型下的截斷半徑dEAM=6。在不改變原子間相對距離的條件下能夠確定3種情況的原子間位置關(guān)系，這里給出其中一組所有原子的位置關(guān)系如圖4和圖5所示。

Figure 4 Stereogram of EAM truncation radius derivation圖4 EAM截斷半徑推導(dǎo)立體圖

Figure 5 Three views of EAM truncation radius derivation圖5 EAM截斷半徑推導(dǎo)三視圖

算法1為空位選取優(yōu)化算法的偽代碼實現(xiàn)。

算法1空位選取優(yōu)化算法

輸入：lattice_constant，vacant，input_i，events，d。

輸出：空位更新后的系統(tǒng)。

1.unitdis←2/lattice_constant；/*晶格常數(shù)轉(zhuǎn)換因子，用于獲取原子坐標(biāo)時將晶格常數(shù)轉(zhuǎn)換為2，從而與截斷半徑的推導(dǎo)過程關(guān)聯(lián)起來*/

2.i←vacant[input_i]；//獲取發(fā)生躍遷事件的空位i

3.proball←sum（events[i].propensity）；/*空位i對應(yīng)的事件列表內(nèi)所有躍遷事件的概率和*/

4.ifproball＞0then

5.threshhold←proball*random（）；

6.j←threshhold；//選取發(fā)生躍遷交換的原子j

7.do_jump（i，j）；//發(fā)生躍遷事件

8.endif

9.d1，d2←0，0；

10.forkinvacantdo

11.round（i，j，k）；/*獲得躍遷原子i，j和查詢空位k的實際坐標(biāo)*/

12.round（i，j，k）*unitdis；/*獲得晶格常數(shù)轉(zhuǎn)換為2后各原子對應(yīng)的坐標(biāo)*/

13.d1←compute_linedis（i，k）；/*計算i與k的距離d1*/

14.d2←compute_linedis（j，k）；/*計算j與k的距離d2*/

15.ifd1≤dord2≤dthen

16.site_propensity（k）；/*更新空位k的躍遷概率*/

17. endif

18.endfor

在403規(guī)模的盒子中進行模擬，分別驗證Pair勢能模型和EAM勢能模型下優(yōu)化算法的正確性。首先用無優(yōu)化的程序，即不考慮截斷半徑，每次躍遷結(jié)束后都對全局空位進行更新的方法對不同空位濃度進行了7組測試，并得到了系統(tǒng)的初始能量和末態(tài)能量。接著用空位選取算法優(yōu)化后的程序在相同的輸入環(huán)境下同樣進行7組測試，并得到了算法優(yōu)化后的系統(tǒng)始末能量，2種模型下的能量對比分別如表2和表3所示。

Table 2 Comparison of system energy before and after Pair potential model optimization表2 Pair勢能模型優(yōu)化前后系統(tǒng)能量對比

Table 3 Comparison of system energy before and after EAM potential model optimization表3 EAM勢能模型優(yōu)化前后系統(tǒng)能量對比

將2種模型無優(yōu)化時的末態(tài)能量畫成曲線的形式，優(yōu)化后的末態(tài)能量畫成點的形式得到圖6，可以看出，優(yōu)化后的點均落在無優(yōu)化時的曲線上面，即2種模型在初始能量相同的情況下，不同空位濃度下算法優(yōu)化后的末態(tài)能量均與無優(yōu)化時的末態(tài)能量保持一致，優(yōu)化后的程序并沒有影響模擬的演化進程，保證了算法優(yōu)化后模擬結(jié)果的正確性。

Figure 6 Final energy comparison before and after optimization of two models圖6 2種模型各自優(yōu)化前后末態(tài)能量對比

接著為了測試優(yōu)化算法的效率提升，本文分別對Pair勢能模型和EAM勢能模型在無優(yōu)化和空位選取優(yōu)化算法下進行了不同空位濃度的測試，并統(tǒng)計了模擬用時。另外，本文選用了動力學(xué)模擬軟件Lakimoca[9]當(dāng)中的截斷半徑數(shù)值（d=10.5）進行空位選取，同樣對2種勢能模型下不同空位濃度的體系進行測試并得到相應(yīng)的模擬用時。最終得到了如圖7和圖8所示的效率提升對比圖。由圖7和圖8可知，隨著空位濃度的增加模擬用時越來越長，EAM勢能模型下的無優(yōu)化模擬用時大幅度多于空位選取優(yōu)化的模擬用時，當(dāng)空位濃度為1.6%時，無優(yōu)化的模擬用時最高達到了空位選取優(yōu)化模擬用時的26.75倍。同時，Lakimoca空位選取方式下的模擬用時相比于無優(yōu)化時有很大的減少，但當(dāng)空位濃度達到1.6%時，仍有超出空位選取優(yōu)化一倍的時間開銷。Pair勢能模型下3種方法的模擬用時與EAM勢能模型下的用時變化具有相同的增大趨勢，但由于Pair勢的計算方式本身比EAM勢更加簡單，因此在相同空位濃度環(huán)境下比EAM勢用時更少。顯然，無論是Pair模型還是EAM模型，空位選取優(yōu)化算法下的模擬都具備更高的效率，尤其當(dāng)空位濃度升高時，能夠避免很大一部分由于空位更新帶來的冗余計算，使得模擬花費的時間開銷進一步減小。

Figure 7 Comparison of Pair potential efficiency improvement圖7 Pair勢效率提升對比

Figure 8 Comparison of EAM potential efficiency improvement圖8 EAM勢效率提升對比

2.3 性能分析

為了得到優(yōu)化后的程序在大規(guī)模并行模擬當(dāng)中表現(xiàn)出的性能數(shù)據(jù)，本文從強可擴展性和弱可擴展性2方面對程序進行了測試與分析。測試環(huán)境的相關(guān)配置參數(shù)如表4所示。

Table 4 Test environment related parameters表4 測試環(huán)境相關(guān)參數(shù)

強可擴展性的測試是指固定模擬體系規(guī)模，只增加CPU核數(shù)，并觀察加速比和并行效率的變化。本文將模擬體系規(guī)模固定為10003，初始化模擬時間為t=1×10-4s，并設(shè)置用來控制自適應(yīng)同步算法的時間間隔參數(shù)nstop=10（nstop參數(shù)為每個子域設(shè)定一個可發(fā)生事件的平均值，結(jié)合所有子域發(fā)生事件的概率均值來計算子域的時間增量）。在溶質(zhì)Cu濃度為1.34%，空位濃度為0.08%，模擬溫度T=573 K的體系環(huán)境中，將CPU核數(shù)從1 000增加到5 000，分別對Pair勢和EAM勢進行算例測試，得到各組算例的模擬用時，并以1 000核算例為基準計算了其他各組算例的加速比和并行效率，測試結(jié)果分別如表5和表6所示。

Table 5 Test results of Pair potential strong scalability（Scale of simulation system：10003）表5 Pair勢強可擴展性測試結(jié)果（模擬體系規(guī)模：10003）

Table 6 Test results of EAM potential strong scalability（Scale of simulation system：10003）表6 EAM勢強可擴展性測試結(jié)果（模擬體系規(guī)模：10003）

2種勢能模型并行效率的變化過程如圖9所示。結(jié)合表5和表6可以看出，當(dāng)CPU核數(shù)從1 000增加到3 000時2種勢能模型的模擬用時下降可觀，呈現(xiàn)出不錯的加速比和并行效率。當(dāng)CPU核數(shù)進一步增加到5 000時加速比增長速度繼續(xù)放緩，其中EAM勢在5 000核時的加速比已經(jīng)小于4 000核時的加速比，并且整體的并行效率也呈現(xiàn)下滑。這是因為CPU核數(shù)增加使得進程間的通信工作量變大，通信時間開銷逐漸占有更高比重，使得并行效率呈現(xiàn)下滑的趨勢。5 000核的測試算例中，Pair勢的并行效率為52%，EAM勢的并行效率為49%，在大規(guī)模體系并行模擬中總體表現(xiàn)出良好的強可擴展性。

Figure 9 Parallel efficiencies of strong scalability of two models圖9 2種模型的強可擴展性的并行效率

弱可擴展性的測試是指保持單個CPU核的工作負載不變，在增大模擬體系規(guī)模的同時以相同的比例增大CPU核數(shù)，并觀察并行效率的變化過程。弱可擴展性的測試采用與強可擴展性測試當(dāng)中相同溶質(zhì)濃度、空位濃度和模擬溫度的體系環(huán)境，將模擬體系規(guī)模從2003逐漸增大到10003，同時保持CPU核數(shù)以相同比例增長，分別對Pair勢和EAM勢進行測試，得到了各組算例的模擬用時，并以2003模擬規(guī)模的算例為基準計算了其他各組算例的并行效率，測試結(jié)果分別如表7和表8所示。

Table 7 Test results of Pair potential weak scalability表7 Pair勢弱可擴展性測試結(jié)果

弱可擴展性測試中2種勢能模型的并行效率隨CPU核數(shù)增加的變化過程如圖10所示。結(jié)合表7和表8可以看出，隨著模擬體系規(guī)模和CPU核數(shù)的同比例增加，模擬用時并沒有保持相對恒定而是逐漸增加，通信時間增加帶來的負面影響變大。在5 000核測試當(dāng)中Pair勢的并行效率為0.60，EAM勢的并行效率為0.68，仍保持著較好的并行性能，因此程序在大規(guī)模體系并行模擬中總體具有良好的弱可擴展性。

Table 8 Test results of EAM potential weak scalability表8 EAM勢弱可擴展性測試結(jié)果

Figure 10 Parallel efficiencies of weak scalability of two models圖10 2種模型的弱可擴展性的并行效率

3 數(shù)值模擬與分析

模擬采用自主開發(fā)的Open KMC軟件，相比于傳統(tǒng)的原子動力學(xué)蒙特卡洛（AKMC）串行方法的效率在大規(guī)模仿真中會受到模擬時間過長和內(nèi)存溢出的限制，Open KMC在AKMC的基礎(chǔ)上將模擬體系分解成多個進程，同時將每個進程的工作負載劃分為更細粒度的扇區(qū)，以此避免進程間的邊界沖突。利用眾核處理器，Open KMC可以實現(xiàn)和加速上億個原子的大規(guī)模KMC并行仿真。本文在文獻[7]的Open KMC程序基礎(chǔ)上實現(xiàn)了2.2節(jié)中的計算優(yōu)化，在不影響模擬結(jié)果的情況下使得模擬速度得到進一步提高，并將優(yōu)化后的Open KMC程序應(yīng)用于接下來的大規(guī)模并行模擬當(dāng)中。另外，文獻[7]采用Open KMC基于Pair勢能模型進行了大規(guī)模并行模擬，本文除了采用Pair勢能模型外，還將加入模擬結(jié)果更為可靠的EAM勢能模型做分析和比較。并行模擬所采用的測試環(huán)境的各項參數(shù)如表4所示。

3.1 并行模擬

首先，構(gòu)造規(guī)模為5003的bcc晶格Fe基體，原子總數(shù)為2.5×108。α-Fe的晶格常數(shù)為2.855?。溶劑原子是Fe，溶質(zhì)原子是占比為1.34%的Cu，體系中引入的空位點缺陷濃度為4×10-4%，模擬溫度T=573 K。在保證內(nèi)存不會溢出的前提下，為了能夠充分利用計算資源，本次模擬采用單節(jié)點64核進行并行模擬。根據(jù)Vincent等人[14]的工作，本文將模擬時間進行了重新標(biāo)度，從而獲得相應(yīng)的物理時間尺度。其時間標(biāo)度轉(zhuǎn)換公式如下所示：

在模擬剛開始階段體系中孤立的Cu原子數(shù)目為3 006 580，體系引入的1 000個空位隨機分布在并行模擬的64個進程當(dāng)中。圖11顯示了Pair勢能模型與EAM勢能模型下體系內(nèi)孤立銅原子的下降過程，模擬對應(yīng)的實際物理時間已經(jīng)通過式（6）和式（7）縮放為蒙特卡洛模擬時間。由于Pair勢能模型下的曲線在模擬初期下降很陡，圖11中右上角的小圖更好地顯示了2種模型在模擬的最早時間內(nèi)孤立銅原子的變化過程。2種模型體系中孤立Cu原子的數(shù)目均在一開始迅速下降，接著下降速度逐漸放緩，最終保持著個位數(shù)的下降趨勢。顯然Pair勢相比EAM勢的模擬結(jié)果仍有差距，但考慮到2種模型下的結(jié)果呈現(xiàn)相同的趨勢，并且由于Pair勢能更小的計算量使得模擬用時大大減少，因此利用Pair勢能模型進行數(shù)值模擬與分析在多數(shù)情況下也是較好的選擇。同時，根據(jù)本次實驗的模擬結(jié)果也可以對Pair勢能模型參數(shù)做進一步的調(diào)整，從而使得模擬結(jié)果更加能接近EAM勢能模型的結(jié)果。

Figure 11 Decline curve of the number of isolated Cu atoms圖11 孤立Cu原子數(shù)量下降曲線

利用EAM勢能模型下的模擬結(jié)果對實驗現(xiàn)象進行分析。經(jīng)統(tǒng)計在整個模擬過程的前0.1 s蒙特卡洛模擬時間里，體系內(nèi)已有超過50%的孤立Cu原子形成了團簇。用時間標(biāo)度公式換算后相當(dāng)于Fe-Cu合金服役工作65年。利用粒子可視化軟件Ovito對0.1 s蒙特卡洛模擬時間下的形核結(jié)果進行輸出。此時體系中已經(jīng)析出大量如圖12和圖13所示的富Cu團簇和Cu-Va的復(fù)合體，其中淺色的圓球代表空位點缺陷，深色的圓球代表Cu原子。

Figure 12 Pure Cu clusters in the system圖12 體系內(nèi)的純Cu團簇

Figure 13 Cu-Va complex in the system圖13 體系內(nèi)的Cu-Va復(fù)合體

為了驗證體系中Cu和空位的富集現(xiàn)象，計算了原子間短程有序SRO（Short Range Order）參數(shù)值[15]，其中有序參數(shù)定義為：

Table 9 SRO value____表9 SRO值

將體系中每幀狀態(tài)下對應(yīng)的最大純Cu團簇尺寸、最大Cu-Va復(fù)合體（空位數(shù)不少于10）尺寸統(tǒng)計輸出如圖14所示。因為體系中空位的數(shù)量遠少于溶質(zhì)Cu原子的數(shù)量，在模擬的前期純Cu團簇更易形成。但是，隨著模擬的進行，空位逐漸聚集并與Cu原子組合形成Cu-Va復(fù)合體團簇，復(fù)合體尺寸反超純Cu團簇尺寸，成為了體系內(nèi)尺寸更大的團簇。這是由于單空位向空位團簇遷移可以使系統(tǒng)能量降低，因此形成更大更穩(wěn)定的團簇；同時，Cu-Va間的結(jié)合能大于Cu-Cu間的結(jié)合能，組成的Cu-Va復(fù)合體可以進一步降低體系的能量，這也是Cu-Va復(fù)合體團簇在體系中可以穩(wěn)定存活的主要原因。

Figure 14 Comparison of maximum cluster size圖14 最大團簇尺寸比較

在溶質(zhì)的擴散聚集過程當(dāng)中綜合了空位形成能和遷移能的作用，以往的AKMC模擬大多在體系中只放入一個空位，并且體系的規(guī)模都很小，例如1×106個原子。而本文模擬結(jié)果表明，當(dāng)體系中空位數(shù)量不止1個時，體系中出現(xiàn)了很多較大的Cu-Va復(fù)合體團簇，這說明傳統(tǒng)的單空位模擬方案未必合理，我們不能忽略空位增加對團簇形成產(chǎn)生的影響。本文模擬的體系中包含1 000個空位，為了研究空位數(shù)量對團簇形成產(chǎn)生的影響，本文在與之前相同的體系環(huán)境中分別引入500，800，1 200，1 600和2 000個空位進行了5組對比測試。首先是體系中孤立Cu原子的數(shù)目變化曲線對比，如圖15所示，可以看出，在沉淀析出過程趨于平穩(wěn)之前空位數(shù)量的增多加快了體系內(nèi)Cu原子的聚集過程。進一步地，為了研究空位數(shù)量對團簇的數(shù)量密度和尺寸產(chǎn)生的影響，提取了蒙特卡洛模擬時間0.1 s下各組測試中最大的團簇，發(fā)現(xiàn)它們的尺寸均不超過40個原子。

Figure 15 Decline curve of isolated Cu atoms in different vacancy number systems圖15 不同空位數(shù)量體系下的孤立Cu原子下降曲線

本文假設(shè)體系中尺寸≥25個原子的團簇為較大團簇，將各組測試環(huán)境中全部團簇（尺寸≥5個原子）的數(shù)量密度和尺寸≥25個原子的團簇的數(shù)量密度統(tǒng)計輸出如圖16和圖17所示。圖16和圖17的縱坐標(biāo)表示單位體積內(nèi)的團簇數(shù)量?？梢钥吹?，隨著空位數(shù)量的增加，體系內(nèi)總團簇的數(shù)量密度相對穩(wěn)定且處于同一個數(shù)量級，沒有較大差異。與此同時空位越多，體系內(nèi)的較大團簇數(shù)量密度有明顯提高。此時體系內(nèi)的細小團簇會出現(xiàn)分解，相對較大的團簇通過吸收細小團簇分解出的Cu原子來保持尺寸增長。而空位增多能夠促進這種團簇粗化現(xiàn)象，進而形成更多較大尺寸的沉淀。

3.2 實驗比較

為了與真實的實驗過程進行比較，本文在403的體系下對663 K～773 K的溫度進行了一系列熱時效模擬。其中溶劑為Fe，溶質(zhì)為占比1.34%的Cu，引入占比8×10-4%的空位濃度。并與Lê等人[5]的實驗結(jié)果進行了比較。本文對模擬的演化過程進行了蒙特卡洛模擬時間10 s的跟蹤，利用式（6）和式（7）轉(zhuǎn)換得到不同溫度下實際演化的物理時間real_time=105s，并分析了這段時間的溶質(zhì)析出過程。105s對應(yīng)的蒙特卡洛模擬時間分別為：溫度為663 K時tMC_Pair=3.2×10-4s和tMC_EAM=9.53×10-3s；溫度為693 K時tMC_Pair=1.3×10-3s和tMC_EAM=3.3×10-2s；溫度為733 K時tMC_Pair=7.04×10-3s和tMC_EAM=0.15 s；溫度為773 K時tMC_Pair=3.2×10-2s和tMC_EAM=0.59 s。按照式（9）計算沉淀推進系數(shù)，其中NX（0）是體系初始化時的孤立銅原子數(shù)量，NX（t）是時間t時體系內(nèi)孤立銅原子的數(shù)量，NX（∞）為無窮遠處的體系內(nèi)孤立銅原子數(shù)量。

Figure 16 Evolution curve of total clusters number density under different vacancy number systems圖16 不同空位數(shù)量體系下的總團簇數(shù)量密度演化曲線

Figure 17 Evolution curve of number density of clusters with large size（≥25）圖17 較大尺寸（≥25）團簇數(shù)量密度演化曲線

ζ（t）=（NX（0）-NX（t））/（NX（0）-NX（∞））

（9）

圖18顯示了4種溫度下的Pair勢能模型和EAM勢能模型分別經(jīng)歷蒙特卡洛模擬時間10 s演化的體系內(nèi)孤立銅原子的下降過程。本文中NX（∞）選用蒙特卡洛模擬時間10 s時的孤立銅原子數(shù)量。將計算后的結(jié)果與電阻率測量得到的真實實驗值進行了比較，結(jié)果如圖19所示。

顯然，在熱時效的模擬過程中，Cu的析出是漸進的。在這一部分工作中，文獻[7]中的模擬結(jié)果與實驗值的對比顯示出了很好的吻合程度。而本文中的模擬結(jié)果與實驗值曲線則有明顯的區(qū)別。這是由于之前的計算并沒有將NX（∞）取為理想的無窮遠處的孤立銅原子數(shù)量，而是在模擬的體系尚未達到相對穩(wěn)定的狀態(tài)時取最后時刻的孤立銅原子數(shù)量作為NX（∞）并代入計算，導(dǎo)致計算出的推進因子并不能實際反映模擬過程中的沉淀推進過程。本文中的沉淀推進過程與實驗值的偏差來自簡化模型的定量偏差，例如系統(tǒng)內(nèi)只有一種類型的缺陷，即空位點缺陷，以及原子間缺乏長程相互作用。但是，模擬曲線至少可以定性地描述體系隨時間增加的演化過程，并且當(dāng)溫度升高時，EAM勢能模型的擬合結(jié)果與實驗值更加接近，在溫度為773 K時，所得模擬結(jié)果與實驗值基本吻合。簇粗化，形成尺寸較大的沉淀。最后，本文計算了4種不同溫度下模擬的沉淀推進結(jié)果，并與實驗值計算出的結(jié)果進行了比較，在一定程度上本文的沉淀推進過程與實驗結(jié)果具有相同的趨勢，并且當(dāng)溫度升高時具備更好的擬合效果。

Figure 18 Evolution of isolated Cu atoms in Pair potential model and EAM potential model at different temperatures圖18 不同溫度下Pair勢能模型與EAM勢能模型的孤立Cu原子演化過程

Figure 19 Evolution of precipitation advance factor in Pair potential model and EAM potential model at different temperatures圖19 不同溫度下Pair勢能模型與EAM勢能模型沉淀推進因子演化過程

4 結(jié)束語

本文基于Pair勢能模型和EAM勢能模型，針對Open KMC軟件中原子躍遷事件發(fā)生后的空位更新問題進行了算法優(yōu)化，在保證模擬結(jié)果與未優(yōu)化時一致的情況下，縮短了模擬用時，并且當(dāng)空位濃度升高時具有更好的效率提升。同時，對優(yōu)化后的程序做了大規(guī)模并行模擬的性能分析，其結(jié)果顯示程序具備良好的強可擴展性和弱可擴展性。進而利用優(yōu)化后的Open KMC軟件并行模擬了熱時效下大規(guī)模體系Fe-Cu合金中引入多空位后的團簇析出過程。模擬結(jié)果表明：（1）本文方法預(yù)測了體系中富Cu沉淀和Cu-Va復(fù)合體的析出，并通過計算SRO參數(shù)驗證了原子聚集現(xiàn)象。（2）統(tǒng)計了體系內(nèi)最大團簇尺寸的變化過程。發(fā)現(xiàn)盡管引入空位點缺陷的數(shù)量遠小于溶質(zhì)Cu原子的數(shù)量，但體系內(nèi)形成的Cu-Va復(fù)合體相比于純Cu團簇更易成為體系尺寸更大的團簇。（3）在相同的體系環(huán)境下，增大引入合金中的空位數(shù)量，可以加快沉淀的析出過程。（4）空位數(shù)量的增加不會對體系內(nèi)總的團簇數(shù)量密度產(chǎn)生明顯的影響，但是能夠促進團