原子動力學蒙特卡洛程序Open KMC在反應堆壓力容器鋼缺陷損傷研究中的優化與應用*

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

（1.太原理工大學大數據學院，山西 太原030024；

2.中國科學院計算技術研究所計算機體系結構國家重點實驗室，北京100190；

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

1 引言

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

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

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

2 程序優化

2.1 基于空位躍遷的KMC方法

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

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

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

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

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

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

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

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

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

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

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

2.2 計算優化

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

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

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

Figure 2 Stereogram of Pair truncation radius derivation圖2 Pair截斷半徑推導立體圖

Figure 3 Three views of Pair truncation radius derivation圖3 Pair截斷半徑推導三視圖

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

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

Figure 4 Stereogram of EAM truncation radius derivation圖4 EAM截斷半徑推導立體圖

Figure 5 Three views of EAM truncation radius derivation圖5 EAM截斷半徑推導三視圖

算法1為空位選取優化算法的偽代碼實現。

算法1空位選取優化算法

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

輸出：空位更新后的系統。

1.unitdis←2/lattice_constant；/*晶格常數轉換因子，用于獲取原子坐標時將晶格常數轉換為2，從而與截斷半徑的推導過程關聯起來*/

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

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

4.ifproball＞0then

5.threshhold←proball*random（）；

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

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

8.endif

9.d1，d2←0，0；

10.forkinvacantdo

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

12.round（i，j，k）*unitdis；/*獲得晶格常數轉換為2后各原子對應的坐標*/

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規模的盒子中進行模擬，分別驗證Pair勢能模型和EAM勢能模型下優化算法的正確性。首先用無優化的程序，即不考慮截斷半徑，每次躍遷結束后都對全局空位進行更新的方法對不同空位濃度進行了7組測試，并得到了系統的初始能量和末態能量。接著用空位選取算法優化后的程序在相同的輸入環境下同樣進行7組測試，并得到了算法優化后的系統始末能量，2種模型下的能量對比分別如表2和表3所示。

Table 2 Comparison of system energy before and after Pair potential model optimization表2 Pair勢能模型優化前后系統能量對比

Table 3 Comparison of system energy before and after EAM potential model optimization表3 EAM勢能模型優化前后系統能量對比

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

Figure 6 Final energy comparison before and after optimization of two models圖6 2種模型各自優化前后末態能量對比

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

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

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

2.3 性能分析

為了得到優化后的程序在大規模并行模擬當中表現出的性能數據，本文從強可擴展性和弱可擴展性2方面對程序進行了測試與分析。測試環境的相關配置參數如表4所示。

Table 4 Test environment related parameters表4 測試環境相關參數

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

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

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

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

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

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

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

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

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

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

3 數值模擬與分析

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

3.1 并行模擬

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

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

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

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

Figure 12 Pure Cu clusters in the system圖12 體系內的純Cu團簇

Figure 13 Cu-Va complex in the system圖13 體系內的Cu-Va復合體

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

Table 9 SRO value____表9 SRO值

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

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

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

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

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

3.2 實驗比較

為了與真實的實驗過程進行比較，本文在403的體系下對663 K～773 K的溫度進行了一系列熱時效模擬。其中溶劑為Fe，溶質為占比1.34%的Cu，引入占比8×10-4%的空位濃度。并與Lê等人[5]的實驗結果進行了比較。本文對模擬的演化過程進行了蒙特卡洛模擬時間10 s的跟蹤，利用式（6）和式（7）轉換得到不同溫度下實際演化的物理時間real_time=105s，并分析了這段時間的溶質析出過程。105s對應的蒙特卡洛模擬時間分別為：溫度為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）計算沉淀推進系數，其中NX（0）是體系初始化時的孤立銅原子數量，NX（t）是時間t時體系內孤立銅原子的數量，NX（∞）為無窮遠處的體系內孤立銅原子數量。

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

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

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

（9）

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

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

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 結束語

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