蒙特卡羅法模擬壓力容器鋼中的團簇

1、蒙特卡羅法模擬壓力容器鋼中的團簇（材料科學與工程學院，上海200072）摘 要：反應堆壓力容器鋼的脆化可能源于團簇。由于在電鏡下很難識別納米團簇，所以通過蒙特卡羅與分子動力學，模擬團簇與位錯的交互作用。模擬結果表明：空位簇阻礙作用最強，其余團簇差異微??；Cu-Ni-v簇的阻力與位錯線的取向有關；小尺寸Cu，Cu-v或Cu-Ni-v簇對位錯阻礙作用較弱； 將Ni替換為Cu原子不會導致強度的任何顯著增加。關鍵詞：壓力容器；團簇；蒙特卡羅模擬；位錯Monte Carlo Simulation of Clusters in Pressure Vessel Steel (School of Materi

2、al Science and Engineering, Shanghai 200072, China)Abstract: The embrittlement of the reactor pressure vessel steel may be due to clusters. Because it is difficult to identify nanoclusters under electron microscopy the interaction between clusters and dislocations is simulated by Monte Carlo and mol

3、ecular dynamics. The results show that the latent cluster has the strongest effect and the other clusters are small. The resistance of Cu-Ni-v cluster is related to the orientation of dislocation line. The Cu，Cu-v or Cu-Ni-v cluster The hindrance is weak; replacing Ni with Cu does not result in any

4、significant increase in strength.Key words：Pressure vessel; cluster; Monte Carlo simulation; dislocation1. 團簇團簇1一般指由幾個甚至上萬個原子或分子形成的相對穩(wěn)定的聚集體，介于孤立原子分子與宏觀物質之間，用無機分子來描述顯得太大，用小塊固體描述又顯得太小，許多性質既不同于單個原子分子，又不同于固體和液體，也不能用兩者性質的簡單線性外延或內插得到。因此，人們把團簇看成是介于原子、分子與宏觀固體物質之間的物質結構的新層次，是各種物質由原子分子向大塊物質轉變的過渡狀態(tài)。團簇結構的形成是由組成團

5、簇的原子數(shù)目和原子間的相互作用勢決定的2。利用勢函數(shù)計算原子團簇的勢能的最小值，來預測團簇的最低勢能時的結構。 2. 第一原理分子動力學模擬團簇2.1第一原理分子動力學原理通過第一原理分子動力學研究團簇1，一般無法嚴謹求解，因而在實際計算中必須采用近似。一個基本的近似就是基于Born-Oppenheimer近似，將電子的運動和原子核的運動分開，以此描述體系電子結構和總能量。在團簇原子核坐標固定的情況下體系總能量的求解常采用兩種方法，密度泛函理論（Density-Functional Theory，DFT）和緊束縛（Tight-Binding，TB）方法。密度泛函理論從第一原理出發(fā)，利用自洽的方

6、法計算出系統(tǒng)的電子結構，可以定量的反映原子間的價鍵強度。但是第一原理密度泛函方法的計算量非常大，不適合計算原子數(shù)上百的大團簇體系。利用第一原理計算團簇穩(wěn)定構型過程中，并沒有充分考慮團簇電子結構的性質。而緊束縛方法通過一套物理參數(shù)來描述原子間的價鍵關系。由這些物理參數(shù)組裝體系電子部分的大哈密頓矩陣，然后在原子核固定的條件下由量子力學的方法求出體系的電子結構，其計算量比密度泛函小 2-3 個量級。根據第一原理（密度泛函和緊束縛）分子動力學，在團簇原子核坐標固定的情況下，計算出體系的電子結構和總能量。再利用 Hellmann-Feynman 定理對電子總能微分就可以求出每個原子核感受到的力。原子核在

7、該力作用下按照經典力學方程運動，從而可以進行結構優(yōu)化、動力學模擬等計算。2.2第一原理分子動力學的優(yōu)點當團簇包含的原子個數(shù)較少時，可結合基于密度泛函的第一原理分子動力學方法進行計算；而對于更大的團簇體系， 為了進一步減小計算量，還可以使用Tight-Binding 方法1。由于密度泛函方法的計算量非常大，不適合計算大團簇體系，因此在計算大團簇體系的時候我們采用基于第一原理的緊束縛方法。由于緊束縛方法的精度取決于所用物理參數(shù)的有效性，只要有合適的參數(shù)，緊束縛方法就能用較小的計算量研究非常大的團簇體系?；诘谝辉淼木o束縛方法和密度泛函方法相比較，不僅保持了第一原理的精度，而且可以以較小的計算量計

8、算較大的團簇體系。緊束縛方法在對原子數(shù)上百的碳納米管的計算中取得了很好的效果。3. 蒙特卡洛法模擬團簇3.1蒙特卡羅法原理蒙特卡洛法又稱隨機抽樣法，基本原理2：對影響可靠度的隨機變量進行大量的隨機抽樣，然后把大量的抽樣數(shù)據代入功能函數(shù)，去確定結構失效概率式。由概率定義知，某事件的概率可以用大量試驗中該事件發(fā)生的頻率來估算，當樣本容量足夠大時，可以認為該事件的發(fā)生頻率即為其概率。因此，可以先對影響其可靠度的隨機變量進行大量的隨機抽樣，然后把這些抽樣值一組一組地代入功能函數(shù)式，確定功能函數(shù)是否失效，最后從中求得函數(shù)的失效概率。3.2 Metropolis蒙特卡羅法3.2.1 算法原理Metropo

9、lis蒙特卡羅算法是一種權重方法（或稱隨機抽樣方法）。在熱平衡情況下，對于用極限頻率分布ens() 表征的特定統(tǒng)計系綜，利用這種方法可生成系統(tǒng)狀態(tài)3。在化學組分恒定的系綜（正則、微正則和等溫等壓系綜）中，其概率分布是系統(tǒng)哈密頓量的函數(shù)；而在化學組分變化的系綜（巨正則系綜）中，概率分布則是化學勢的函數(shù)。對于由N個經典粒子組成的系統(tǒng)，其熱力學狀態(tài)可以用位置矢量（r1 ，rN ）的三個分量和動量矢量（p1 ，pN ）的三個分量描述。這 6N個分量可以看作矢量在6N維相空間中的分量。3.2.2.正則和微正則系綜的Metropolis方法3 在正則和微正則系綜中，蒙特卡羅方法通過任意地或系統(tǒng)地改變粒子位

10、置可以給出新的組態(tài)。例如，把原子的位置由i變到j就代表系統(tǒng)狀態(tài)的變化，并記為i j，在相空間它表示系統(tǒng)由i點運動到j 點。根據所考慮粒子之間的相互作用，這一點位移將使系統(tǒng)能量由 H( i )變?yōu)镠( j) 。然而，原子的這一位移并不影響系統(tǒng)的組成。計算出新組態(tài)相關的系綜哈密頓量變化H(ij)，即H(ij)= H( j) - H( i ) 式中，H( j)為試探態(tài)的哈密頓量；H( i )為初始態(tài)的哈密頓量。然而，根據下面介紹的定則，可以計算出新的原子排列的情況。如果新的能量值比變化前的能量值小，則該位移使得系統(tǒng)處于較低能量狀態(tài)。這時，變更即被接受且其位移原子保留在新的位置上。如果新的能量值大于前

11、一能量值，則該變更僅以一定的概率ens(j) /ens(i) 被接受，在正則系統(tǒng)中這個概率由玻耳茲曼因子給出，即：ensjensi=exp-Hij式中， =1/kB T。當H(ij) 取正值時，接受組態(tài)變化的概率為：Pij exp-Hij 根據Metropolis等人的理論，若假定在0到1之間由一個隨機數(shù) ，則新的組態(tài)可按照下列規(guī)定確定：當exp- H(ij)時，接受變更，否則不接受變更。如果新的組態(tài)被駁回，亦即沒有被接受，則把原位置記為新位置并記數(shù)一次，然后采用其它隨機選定的原子重復上述過程。3.2.3巨正則系綜的 Metropolis方法3 巨正則系綜中，其基本組分的濃度是沒有限定的。這樣

12、，就可任意選擇一個原子并通過改變一次原子種類而得到一個新的組態(tài)j。顯然，這一過程將影響到系統(tǒng)的化學組成。同樣，組態(tài)的變更將以一定的概率被接受。在巨正則系綜，因其組成變化而引起的能量變化U(ij)，也可由計算得到。對新的組態(tài)可按照下列標準進行判斷：如果得到的新的能量小于原來的能量，則組態(tài)的這一變更被接受；如果能量大于原來的能量，則被一定的概率接受。當U(ij)為正值時，被接受的概率：Pij exp-Uij （2.16）式中，U 表示混合能及化學勢該變量之和。判斷的依據與正則和微正則的方式相同。3.3蒙特卡羅法優(yōu)點MC方法，主要適用于原子尺度和顯微尺度的模擬。MC 法用隨機數(shù)來控制粒子運動，并使其

13、符合 Boltzmann 分布。MC方法并非只能用來解決包含隨機過程的問題，例如:對計算定積分的問題可將其轉換成相關的隨機過程， 然后用Monte Carlo方法進行求解。MC與傳統(tǒng)數(shù)學方法相比，具有直觀性強，簡便易行的優(yōu)點，能處理一些其他方法無法解決的問題，并且容易在計算機上實現(xiàn)，在很大程度上代替許多大型、難以實現(xiàn)的復雜實驗。3.4蒙特卡羅法在材料中的應用Monte Carlo方法在材料中的應用主要分為4：經典例子體系中的應用，如團簇結構、薄膜生長、合金偏析和高分子材料等；經典自旋體系，如Ising自旋玻璃、一維下的熱平衡、二維下的二級相變。圖1 MC法在經典例子體系中的應用3.5 蒙特卡羅

14、模擬的步驟： 根據欲研究的物理系統(tǒng)的性質，建立能夠描述該系統(tǒng)特性的理論模型，導出該模型的某些特征量的概率密度函數(shù)（即構造或描述問題的概率過程）； 從概率密度函數(shù)出發(fā)進行隨機抽樣，實現(xiàn)從已知概率分布的抽樣，得到特征量的一些模擬結果；有了明確的概率過程后，為了實現(xiàn)過程的數(shù)值模擬，必須實現(xiàn)從已知概率分布的隨機數(shù)的抽樣，進行大量的隨機模擬實驗，從中獲得隨機變量的大量試驗值。產生已知概率分布的隨機變量，是實現(xiàn)MC方法的關鍵步驟，其中最基本的是（0，1）均勻分布。 對模擬結果進行分析總結，預言物理系統(tǒng)的某些特性。 模擬結果的檢驗。對于MC方法模擬納米團簇的具體算法4見圖2。圖2 MC法的具體算法流程圖4.

15、 具體應用：刃型位錯與bcc鐵中的Cu-Ni空位簇的相互作用54.1案例背景反應堆壓力容器鋼脆化的硬化原因6：銅高度不溶于鐵，在中子照射下導致團簇（2nm）形成，富銅的析出物引起輻照硬化。納米空位簇和位錯環(huán)造成的基質損傷。由于不能檢測低于TEM分辨率的位錯環(huán)，無法論證不同類型鋼中的位錯環(huán)的存在7，位錯與特定合金元素的關聯(lián)不能通過實驗明確地建立。所以通過蒙特卡羅與分子動力學，考慮它們與移動位錯的相互作用。研究刃型位錯與包含銅原子，鎳原子和空位的納米團簇以及純納米空位簇和純銅簇的相互作用6。bcc鐵中的空位簇、Cu簇與位錯的相互作用MD研究，已經表明：只有大的純Cu團簇（ 3nm）作為位錯的強障礙

16、物；可見的位錯環(huán)（2nm）可能作為位錯的強障礙；低于TEM分辨率的位錯環(huán)在與位錯相互作用時自發(fā)地改變它們的伯氏矢量并且被立即吸收；溫度對位錯與1nm以下團簇的相互作用幾乎沒有影響；空位簇是比純Cu團簇更強的障礙物，小于1nm的空位簇容易被位錯剪切。4.2模擬方法為了確定給定組成的每個簇中的溶質和空位的能量最有利的排列，我們使用Lattice Metropolis Monte Carlo（LMMC）采樣技術獲得復雜簇的原子構型。LMMC計算在包含根據bcc結構排列的2000個原子的立方體盒中進行，其中主軸沿著100，010和001方向，并且空位，Ni和Cu原子初始隨機分布。不同晶格位置處不同類型

17、的原子交換導致總能量的變化，在晶體中產生的三元原子間勢。模擬中，拒絕總能量增加的步驟（即原子位置交換），運行總能量減少的步驟，直到晶體的總能量收斂到恒定值。隨后，將團簇模型輸出到晶體模擬盒中的位錯滑移面上?；谔剿鞑煌娜苜|/空位比率以及用Ni原子替換部分Cu原子時的化學效應的想法，選擇了許多團簇模型，見表1，引入純Cu簇和空位簇，以便參考比較。用于確定簇中的原子和空位的空間分布的程序如下：使用原子離析分析，耦合結構最近鄰分析，以及通過選擇具有高勢能的原子來監(jiān)測模擬期間位錯核、團簇的位置和結構。表1 Cu-Ni-v簇的構型及其總結合能序號空位數(shù)銅原子鎳原子團簇大小能量差11020030-13.

18、4421040050-20.5331060070-27.1340151530-6.012510301050-18.59610303070-18.807700070-73.368500050-46.679300030-23.9210070070-21.4611050050-13.8212030030-7.5晶體(bcc)模擬盒的主軸x，y和z分別沿著111、2和10?；泼鏋?110)，位錯線與y方向平行，并且具有平行于x軸的伯氏矢量=12 111。 沿x和y方向應用周期邊界條件。沿著z，盒子被分成三部分：上部和下部由幾個原子平面組成，其中原子剛性地固定在其原始位置，而內部區(qū)域中的原子在MD運行

19、期間自由移動。 位錯的動力通過在x方向上的塊體剛性移動而產生8。變形引起的剪切應力為s= Fx/Axy，其中Fx是由于內部區(qū)域中的原子在下部剛性塊上沿x方向的總力，Axy是x-y截面面積。為了為該彎曲提供足夠的空間，MD箱的內部區(qū)域的尺寸為沿著x，y和z有1003、306和252個原子平面。因此，盒體積為約252110 nm3，并含有約50萬個移動原子。首先放置直徑為約1nm的Cu-Ni-v或Cu-v簇，其中心位于位錯滑移平面上，離位錯線約5nm，如圖3所示。圖3 MD模擬設置的示意圖；刃型位錯位于含有約50萬個原子的晶體中用MC獲得的Cu-Ni-v簇放在滑移面上；團簇包含70個Cu原子和10

20、個空位4.3結果分析4.3.1 互動機制和一般趨勢最初，位錯滑向團簇，產生交互作用，應力增加，位錯線彎曲，最終從團簇斷開，在團簇出口表面上留下臺階6。圖4可以看出 ，空位簇阻礙作用最強，其余團簇差異微小。觀察到位錯和沉淀物之間的反應通過相同的機理發(fā)生，而與團簇的化學組成或尺寸無關。圖4 刃型位錯通過空位簇，Cu70，Cu60v10和Cu30Ni30v10簇的剪切應力圖5中顯示了純納米空隙和純Cu沉淀的斷裂處位錯的關鍵構型。圖b，與空位簇的相互作用導致位錯偶極子的形成，并且臨界角接近零；Cu團簇（cos77= 0.21） 空位簇（cos18= 0.95）。在與空位簇的相互作用中，觀察到從團簇延出

21、的兩個螺旋臂爬升，因此位錯運動中吸收了一些空位?？瘴坏目倲?shù)與納米空隙直徑成正比;納米空隙越大，空位被吸收的越多。圖5 位錯與V70和Cu70簇在（110）平面交互后的位錯核（a）顯示了模擬框的實際尺寸；（b）障礙物附近的位錯的構造被放大4.3.2 Cu-Ni-v簇的阻力發(fā)現(xiàn)剪切給定的Cu-Ni-v簇所需的臨界應力取決于位錯線取向。 圖6中提供了取向對克服Cu60v10簇所需的臨界應力的影響的實例。即使臨界應力非常小，我們可以看出，一個特定的方向提供比其他兩個明顯更高的阻力。然而，該差異傾向于隨著溫度的升高而消失。即使在與Cu-Ni-v簇的反應中沒有觀察到螺旋偶極子，存在于簇中的部分空位被刃型位

22、錯吸收。 發(fā)現(xiàn)實際吸附的空位數(shù)取決于簇的具體取向，其大小和模擬溫度。 通常，溫度越高，吸附的空位越多。 因此，這些團簇的有效尺寸在與刃型位錯的相互作用之后減小，團簇發(fā)生原子重排。通過提高溫度，Cu-Ni-v簇的臨界應力也降低。然而，由于由熱原子振動引起的分辨應力的波動，這些團簇的缺點導致這種減少難以發(fā)現(xiàn)。 圖6，取向對克服Cu60v10簇所需的臨界應力的影響。即使臨界應力非常小，11 表現(xiàn)出更高的阻力。然而，該差異傾向于隨著溫度的升高而消失。圖6 位錯通過Cu60v10簇的剪切應力；三個曲線對應于微晶旋轉后新的x軸分別與100,010或001方向一致時的情況圖7a和b分別顯示了在此研究的最小和最大團簇的數(shù)據。 在所有溫度下，所有Cu-Ni-v簇的強度小于或等于相同尺寸的純Cu簇的強度。 因此可以得出結論：所有小的Cu，或Cu-v或Cu-Ni-v簇對位錯阻礙作用較弱， 將Ni替換為Cu原子不會導致強度的任何顯著增加。圖7 模擬中不同簇對應的剪切應力5. 結論使用LMMC技術獲得低于電子顯微鏡