鋯級聯(lián)的分子動力學模擬.

1、鋯輻照級聯(lián)的分子動力學模擬研究小組成員：韓維強（主持人）、司宇、張錦文指導老師：侯懷宇 2012.32013.91 緒論1.1研究背景和意義能源,信息,材料,通常被稱為新科技革命的三大支柱。人類文明的進程和能源的發(fā)展革命息息相關。自工業(yè)革命以來，以化石燃料作為主要能源的時代已經(jīng)持續(xù)了兩百多年。當今世界，人類對能源的需求與日俱增，而地球上探明的可使用的化石燃料如煤炭，石油，天然氣和其他可燃氣體迅速減少，因而尋找一種可替代化石能源的其他能源是迫切而必然的選擇。半個多世紀以來，核電，火電和風力發(fā)電一直作為當今世界電力的三大支柱發(fā)揮著不可代替的作用。核能，以其高效，清潔，安全的優(yōu)點備受青睞。在我國，目

2、前已建成秦山核電站、廣東大亞核電站、江蘇田灣核電站和嶺澳核電站四座核電站，為我國電力事業(yè)發(fā)揮巨大作用1。核材料應具備以下條件：a. 符合核條件。材料的吸收中子截面要??；雜質(zhì)含量滿足核級純；低活性，感生放射性小，剩余發(fā)熱量小，半衰期短。b. 運行性能。應具備良好的物理性能；高抗輻射性能；與冷卻劑相容性好。c. 經(jīng)濟適用性。材料具備工業(yè)生產(chǎn)能力，可批量生產(chǎn)；低成本；焊接性好?；诟黜椥阅芸紤]，工業(yè)應用的核材料主要有鋁及鋁合金、鎂合金、不銹鋼、鋯合金。鋯合金比不銹鋼的熔點高300-400,熱脹系數(shù)小2/3,熱導率高18%,熱中子吸收截面小一個量級;力學性能、加工性能好,同UO2相容性好,尤其是對30

3、0-400的高溫水、高溫蒸汽也具有良好的抗腐蝕性能和足夠的熱強性?；诖?鋯被選作反應堆燃料。元件的包殼材料和結構材料,如壓力管、容器管和定位格架等,是“核安全最堅韌防線”。核反應堆的必要條件除了等離子物理外，最為重要的就是壁材料。所以說，核反應堆的誕生基于鋯的發(fā)展、成熟2。1.2國內(nèi)外輻照級聯(lián)模擬研究現(xiàn)狀目前，輻照損傷的模擬研究方法主要有兩種：實驗模擬和計算機模擬。在實驗研究方面，高能電子輻照3、中子輻照4、離子輻照5-7和等離子輻照8等，是國內(nèi)外模擬材料在核反應堆內(nèi)的中子輻照，對核材料結構和性質(zhì)的探究，欲發(fā)現(xiàn)其規(guī)律時常用的輻照模擬技術。上述模擬方法各有特點，如離子輻照具有實驗條件可控等獨特

4、優(yōu)勢、電子輻照不會引入雜質(zhì)原子,因而可以利用它們的優(yōu)點作綜合分析研究。在核材料的研究領域中,計算機模擬技術得到廣泛的應用。相對實驗方法來說,計算機模擬方法實現(xiàn)起來比較容易,實驗費用較低,且能夠獲得輻照條件下材料變化的微觀細節(jié),得到許多宏觀實驗中無法獲得的信息,輔以相應的輻照模擬實驗及理論分析,就能對輻照損傷進行更為深入細致的了解。輻照損傷的計算機模擬的主要電對象是核反應材料（如SiC，不銹鋼）和航天器工程材料（如石英玻璃）。1.3本工作研究內(nèi)容本文針對核結構材料的輻照損傷過程,以金屬Zr為研究對象,采用分子動力學(MD)方法對輻照損傷初始過程級聯(lián)碰撞進行計算機模擬,展現(xiàn)了輻照損傷在原子尺度的演

5、變過程,并對此過程的溫度影響因素進行了研究。本文利用分子動力學方法，模擬計算了原子尺度下密排六方結構金屬Zr(組成核結構材料Zr合金的主要元素)的輻照損傷過程,對其缺陷的產(chǎn)生演變機制進行了研究,重點討論了PKA入射能量對材料輻照損傷過程的影響。在一定溫度下，本論文分別模擬了具有不同能量的初級離位原子PKA（1Kev,5Kev,10Kev）入射時的輻照過程。模擬發(fā)現(xiàn)PKA能量增大時,輻照引起的級聯(lián)碰撞過程更為劇烈,這說明更多原子得到了超過離位閉值Ed的能量,從而克服束縛力并從晶體點陣移出。2 分子動力學（MD）模擬方法 2.1 分子動力學模擬的概念分子動力學模擬(moleculardynamic

6、s simulation , MD)是在評估和預測材料結構和性質(zhì)方面模擬原子和分子的一種物質(zhì)微觀領域的重要模擬方法,通過計算機對原子核和電子所構成的多體體系中的微觀粒子之間相互作用和運動進行模擬,在此期間把每一原子核視為在全部其他的原子核和電子所構成的經(jīng)驗勢場的作用下按照牛頓定律進行運動,進而得到體系中粒子的運動軌跡,再按照統(tǒng)計物理的方法計算得出物質(zhì)的結構和性質(zhì)等宏觀性能.。簡而言之即是應用力場及根據(jù)牛頓運動力學原理所發(fā)展的一種計算機模擬方法。經(jīng)典的分子動力學方法是由Alder9和wainwright10在1957年和1959年提出并應用于理想“硬球”模型的,他們發(fā)現(xiàn)了早在1939年就由Kir

7、kwood根據(jù)統(tǒng)計力學預言的“剛性球組成的集合系統(tǒng)會發(fā)生由其液相到結晶相的相轉變”，后來人們稱這種相變?yōu)锳lder相變,此舉開創(chuàng)了用分子動力學模擬方法研究物質(zhì)宏觀性質(zhì)的先例。用分子動力學方法，有三點決定了其結果的可靠性：勢函數(shù)的優(yōu)劣邊界條件的合力配置以及計算機的運算速度。勢函數(shù)描述了原子間的相互作用，邊界條件反映了模擬對象所處的真實環(huán)境，還有高性能的計算機是實現(xiàn)分子動力學模擬的必要保證。2.2 分子動力學方法的原理分子動力學將連續(xù)介質(zhì)看成是由N個原子或分子組成的粒子系統(tǒng),體系中每個粒子都遵循牛頓第二運動定律。MD模擬中體系隨著時間的演變是通過對粒子的運動方程積分實現(xiàn)的"通過數(shù)值求解得

8、到粒子在相空間的運動軌跡,然后由統(tǒng)計物理學原理得出該系統(tǒng)相應的宏觀動態(tài)和靜態(tài)特性。MD模擬是用來計算經(jīng)典多體體系的平衡和傳遞性質(zhì)的一種確定性方法。經(jīng)典是指體系中的粒子的運動都遵從經(jīng)典力學定律。簡單的說,就是MD中處理的多體體系的粒子的運動遵從牛頓方程，即 Fi(t)=miai(t) (2.1)式中，F(xiàn)i(t)是粒子i所受的力，mi是粒子的質(zhì)量，ai(t)是原子i的加速度。原子i所受的力Fi(t)可以直接用勢函數(shù)對坐標ri的一階導數(shù)，即Fi(t)=-dU/dri,其中U為勢函數(shù)。不同于蒙特卡洛方法,分子動力學是確定性的方法,即,原子的初始位置和速度給出后,體系內(nèi)任意時刻的坐標及速度都可以通過數(shù)值

9、求解運動方程來確定。D整個運行過程中的坐標和速度稱為軌跡(trajeetory),數(shù)值解普通微分方程的標準方法為有限差分法。2.3周期性邊界條件及作用勢2.3.1 周期性邊界條件采用分子動力學對材料進行計算機模擬，必須首先對被計算的粒子系綜給定適當邊界條件。這些條件一般分為四種：1）自由邊界條件。常用于大型的自由分子的模擬。2）固定邊界。在所有要計算到的粒子晶胞之外包裹幾層結構相同的、位置固定的粒子,包層厚度大于粒子間相互作用力程范圍。包層部分代表了與運動粒子起作用的宏觀晶體的那一部分。這種邊界條件常用于對點缺陷等性質(zhì)的研究。3）柔性邊界。它允許邊界上粒子有微小移動以反映內(nèi)層粒子對其施加作用力

10、的情況。常用于模擬缺陷延伸情況。4）周期性邊界條件模擬較大系統(tǒng)時,為了消除尺寸效應,常采用此類邊界條件。它使得利用較少粒子來模擬各種系綜成為可能。2.3.2 原子間作用勢在MD模擬中，作用勢函數(shù)的選取是個重頭戲，它對模擬起著決定性的作用。勢函數(shù)是描述原子間作用的函數(shù)，原子間的力場相互作用控制著原子的行為，根本上是決定了材料的性質(zhì)。通俗的說，勢函數(shù)是這樣一個模型：作用在粒子上的力F隨著原子位置的變化而變化。一般而言，主要有對勢、多體勢。對勢是指僅有兩個原子的坐標決定的相互作用，這是僅考慮到原子兩兩作用，而忽略其他原子對它們的作用，半導體，金屬以外其他的無機化合物的相互作用可以用這樣的勢函數(shù)較充分

11、地描述。多體勢的研究對象通常是相互作用較為復雜的多粒子系綜，也是實際研究的大多數(shù)情況。2.4 牛頓方程求解算法 我們需要用積分方法求解，這種情況下，牛頓運動方程可用有限差分方法來求解。有限差分方法的思想是將積分分成很多小步，每一小步固定時間是t，在時間t時刻,作用在每個粒子的力的總和等于它與其它所有粒子的相互作用力的矢量和。根據(jù)此力,我們可以得到此粒子的加速度,結合t時刻的位置與速度,可以得到t+t時刻的位置和速度。這力在此時間間隔期間假定為常數(shù)。作用在新位置上的粒子的力可以求出,然后可以導出t+2t時刻的位置與速度等,以此類似。MD模擬中,一個好的積分算法需要具有以下條件:1,計算速度快;2

12、,需要較小的計算機內(nèi)存;3,允許使用較長的時間步長;4,表現(xiàn)出較好的能量守恒。計算中會涉及到截斷誤差(truncation error)和計算機的舍入誤差(round-of error).前者與選取的積分算法以及時間步長有關,后者與計算機如何實現(xiàn)這個算法以及存儲的精度有關.通常,可以通過減小步長t來減少這兩個誤差。在模擬中,應該根據(jù)問題所需要的精度來選取合適的時間步長。常見的算法有Verlet算法，跳蛙算法，Gear預測-校正法等等。積分步長t較大時，用Verlet算法比較合適，積分步長t較小時，使用預估-校正法誤差比較小。本文使用的就是Gear預測-校正算法。接下來，介紹一下這兩種算法。2.

13、4.1 Verlet 算法 Verlet算法于1967年被提出"在分子動力學中,積分運動方程運用最廣泛的方法是Verlet算法。這種算法運用t時刻的位置和速度及t-t時刻的位置,計算出t+t時刻的位置r(t+t)。 verlet算法的推導可以通過對原子位置r(t)進行Taylor展開到三階: (2.1) (2.2)式中,v是速度,a是加速度,b是位置r對時間t的三階導數(shù)。兩式相加,則得出 (2.3)這是verlet算法最基本的形式。式(2.2)是差分公式,其中(t)的奇次項全部被抵消,同時忽略(t4)及其高次項,差分誤差為O(t4)的量級。在模擬中只需要提供原子當前時刻t及前一時刻t

14、-t的位置(當前時刻下a= F/m)就可以得到下一時刻t+t的位置.速度并未出現(xiàn)在verlet算法中,一個簡單的方法是將2.1及2.2兩式相減,即用t+t時刻與t-t時刻的位置差除以2t,即可得到速度的計算公式 (2.4)速度的誤差與時間步長t的二次項相關,誤差在O(t3)的量級,而位置的誤差僅與四次項相關。Verlet算法雖然簡單且存儲要求適度,但有一些缺點:計算過程容易造成精度損失;方程中無顯式速度項,則在下一步位置未確定前,速度項難以求出;非自啟動算法,新位置需由t和t-t的位置得到。針對Verlet算法的這些缺點,人們對Verlet算法做了一些改進,提出了新的算法,例如對速度的處理更好

15、的辦法還有速度Verlet算法和蛙跳算法。2.4.2 Gear預估-校正算法Gear預估-校正是C.W.Gear在1971年提出的,是一類算法的總稱。其核心思想是通過當前時刻的位置及位置對時間的前幾階導數(shù)（一階為速度，二階為加速度）開始，推測下一個時刻的位置及速度。通過預估位置得到的預估加速度與直接預估的加速度之間存在一個偏差；然后通過加速度的差值來校正位置和速度等。在校正的步驟里,需要事先確定對應好不同量的校正參數(shù),因此即使對于同階的算法,不同參數(shù)的選取會產(chǎn)生不同的具體算法。這種方法可分為三步： （2.5）式中,v是速度(位置對時間的一階導數(shù)), a是加速度(位置對時間的二階導數(shù)),b是坐標

16、對時間的三階導數(shù)等。第二步根據(jù)新預測的位置rip(t+t),計算t+t時刻的力F(t+t),然后計算加速度aic（t+t）。這加速度再與由Taylor級數(shù)展開式預測的加速度aip(t+t)進行對比,根據(jù)兩者之差在校正步里用來校正位置與速度項。通過這種校正方法,可以估計預測的加速度誤差為 (2.6)假定預測的量與校正后相差很小,我們可以說,它們互相成正比,這樣校正后的量為 (2.7)Gear確定了一系列的系數(shù)c0,c1,.,展開式在三階微分b(t)后被截斷。采用的系數(shù)的近似值為c0=1/6,cl=5/6,c2=1和c3=1/3。Gear的預估一校正算法需要的存儲量為3(O+l)N,O是應用的最高

17、階微分數(shù),N是原子數(shù)目。2.5 MD模擬的粒子系綜系綜（Ensemble）即系統(tǒng)的綜合，它是統(tǒng)計力學中的一個概念，是由1901年由吉布斯創(chuàng)立完成的。在進行分子動力學模擬中，需要根據(jù)實際條件選定一定的系綜。分子動力學中運用系綜，是為了研究系統(tǒng)的微觀狀態(tài)與宏觀性質(zhì)對應的規(guī)律。分子動力學的研究對象是多粒子系綜，其模擬的粒子數(shù)目有限，而且受到計算機性能的限制，但統(tǒng)計物理規(guī)律仍成立。物質(zhì)世界是由聚集在一起的大量粒子構成,而宏觀物質(zhì)的行為是由微觀粒子間的相互作用決定的。為了研究宏觀體系的演化和結構,基本方法就是構建一個粒子系綜。系綜是一個巨大的系統(tǒng),它是指在一定的宏觀條件下,大量組成、性質(zhì)、尺寸和形狀完全

18、相同的、處于各種運動狀態(tài)的、各自獨立的系統(tǒng)的集合。其中的每個系統(tǒng)各處在某一微觀運動狀態(tài),而且各自獨立。微觀運動狀態(tài)在相空間中構成一個連續(xù)區(qū)域,與微觀量相對應的宏觀量是在一定的宏觀條件下所有可能的運動狀態(tài)的平均值。MD模擬包括平衡態(tài)和非平衡態(tài)模擬。主要系綜包括：1）能量和粒子數(shù)都固定的系綜（N,V,E）：體系與外界不能交換能量，系統(tǒng)原子數(shù)N，體積V，能量E都保持恒定。一般而言，給定能量的精確初始條件無法知曉，為了把系統(tǒng)調(diào)節(jié)到給定的能量，需要先給出一個合理的初始條件，然后對能量進行調(diào)整，使系統(tǒng)達到所要狀態(tài)。一般是對速度V進行特別標度來調(diào)整能量，但這樣會使系統(tǒng)偏離平衡，因此要給足時間來重新讓系統(tǒng)建立

19、平衡。2）原子數(shù)N，體積V，和溫度T都保持不變，總能量為零的系統(tǒng)（N,V,T）：讓系統(tǒng)與虛擬熱裕保持熱平衡狀態(tài)來控制溫度不變。由于溫度和動能有直接關系，所以要固定系統(tǒng)的動能。3）原子數(shù)N，溫度T，壓力P都恒定的等溫等壓系綜：溫度的恒定是通過調(diào)節(jié)系統(tǒng)速度或加一約束力來實現(xiàn)的。P，V是相關量，是通過標定體積V來調(diào)節(jié)壓力值的。2.6 MD運行基本步驟 雖然分子動力學核心思想比較簡單,但實際上分子動力學的實現(xiàn)在很多方面都面臨挑戰(zhàn),如模型的選擇、初始條件的設定,為保證可靠性而進行的各種模擬方案,積分算法的選擇,需要考慮到運動軌跡對初始條件及其他參數(shù)選擇的敏感性,滿足大計算量的要求,圖形顯示與數(shù)據(jù)分析等等

20、11。分子動力學的模擬流程大致如下：1.選取勢函數(shù)，即是建立一個計算模型。勢函數(shù)與物質(zhì)描述密切相關。當勢函數(shù)模型確定時，可以根據(jù)物理學規(guī)律得到模擬中的守恒關系。2.初始化坐標和速度。MD過程中對運動方程的求解，必須知道每個粒子初始位置和速度。（不同的積分方法要求不同的初始條件。）3.趨于平衡的計算。MD模擬中需要設計一個使系統(tǒng)達到一個穩(wěn)定的平衡狀態(tài)的過程。當系統(tǒng)的動能和勢能的總能量在某個值附近波動時，認為已經(jīng)達到平衡。這段達到平衡狀態(tài)的時間成為弛豫時間。4.最后針對模擬對象進行物理量的計算。參與計算的物理量是在模擬的最后進行的。它是沿相空間軌跡求平均得到的。 2.7本文MD模擬采用的軟件情況

21、本文分子動力學模擬采用的軟件是 Lammps。Lammps是分子動力學模擬常用的一款軟件，其源代碼公開，免費下載，可以根據(jù)自己的需要修改lammps代碼，重新編譯；lammps需要在linux環(huán)境下運行，可以串行和并行運算，高移植性的C+語言編寫；lammps是美國能源部下屬的圣地亞國家實驗室發(fā)布的，主要作者是Steve Plimpton。2.7.1 Lammps軟件優(yōu)點可以在Linux并行環(huán)境中運行，體現(xiàn)出強大的計算優(yōu)勢；能模擬上百萬的原子體系，氣態(tài)、液態(tài)或者固態(tài)，在各種系綜下；提供了各種勢函數(shù)可供選擇（只是提供了各個函數(shù)的表達式，具體的參數(shù)是自己找的）。2.7.2 Lammps 運行簡要步

22、驟從文件讀取，可包括拓撲信息構造晶胞復制命令設置力場參數(shù)模擬參數(shù)輸出文件。3 鋯輻照損傷的MD模擬分子動力學的方法雖然能很好地模擬在極端情況下的輻照現(xiàn)象，但是受到很多方面的影響。如果進行一定程度的修改，會使得結果更貼近真實情況。3.1模擬方法3.1.1 本文所用勢函數(shù) 本文采用Mendelev在2007年提出的Zr的EAM勢函數(shù)12。使用的是文獻13中的2號勢函數(shù)。EAM勢函數(shù)的總能量：， （3.1）其下標i，j表示N原子系統(tǒng)中的每個原子， 是成對可能性，是嵌入能量函數(shù)。其中， （3.2）是另外一個成對可能：“密度函數(shù)”為建立一個理想的EAM勢函數(shù)需要給尋找最佳函數(shù)。 本函數(shù)正確描述了HCP到

23、bcc相的過度和液體結構數(shù)據(jù)。作為第一個測試，其可能性被應用在研究HCP和bcc的自擴散，發(fā)現(xiàn)了異常。伴隨自我間遷移的高擴散性與bcc內(nèi)在結構缺陷的形成有關。 3.1.2可變時間步長輻照條件下，原子間的距離會變得非常小（1À甚至更?。?，從而原子間的作用力會非常大，標準的勢函數(shù)很難描述原子間的作用勢，單一的勢函數(shù)會使得誤差增大。因而需要設置這樣一種機制，即當原子受到高能量作用時，時間步長自動縮小，當原子受到低能量作用時，時間步長會自動增大。這樣的機制既保證了準確，又最大限度地控制了模擬的運行時間。在輻照模擬中，常用的算法是依據(jù)計算出的時間步長，對原子在每一個MD時間步長內(nèi)移動的距離加一

24、個限制。一般會給定一個時間步長（常用1fs,0.5fs）在模擬中，我們采用可變時間步長的方法，步長不斷調(diào)整。設定一個最小時間步長，最大時間步長為。每一個時間步長內(nèi)移動的距離限制在一定范圍，最大移動范圍設定為，使得所有原子基于當前原子的速度和力，移動的距離都不超過最大距離。當從一個原子團簇開始模擬或運行一個或更多原子撞擊固體引起級聯(lián)損傷反應的模擬時可能都是有用的。每步之后，下一步的時間步長dt按下列方式計算：對于每個原子，時間步長計算，將導致它的位移在接下來的集成步驟，為其當前的速度和力的作用。由于執(zhí)行此精密計算需要解決一個四次方程，產(chǎn)生一個更簡單的估測。估測是保守的，盡管它可能較小，但原子的位

25、移保證不會超過。3.1.3 缺陷辨別方法及確定空位和離位原子122 A B 圖3-1 完整晶體結構與有缺陷的晶體結構 在輻照模擬分析之前，需要說明缺陷的判別方法。缺陷判別方法，就是對比完整晶體結構和有缺陷的晶體結構，根據(jù)兩結構的不同來尋找缺陷。如圖3-1，完整晶體結構為A，有缺陷的晶體結構為B。B中白圓圈表示空位，右下角的黑圓圈表示間隙原子。在B中的某一原子如果不能在A中找到對應的原子，則認為它是間隙原子。對格點進行循環(huán)檢測：1.每個格子點，確定其所在晶胞的坐標（x，y，z）；2.確定其上下左右的9個晶胞，對每個晶胞內(nèi)的原子進行原子-格點的距離計算；根據(jù)，確定格點的占據(jù)原子（可以多于一個）并確

26、定離格點最近的原子，3.如果離格點最近的原子距離大于，則該點是空位；4.格點占據(jù)者中，并非最近原子的，就是與該點相聯(lián)系的間隙原子。上法已確定一部分間隙原子，若原子并非任何一個格點的占據(jù)者，那它也是間隙原子。3.1.4 PKA速度的確定本文使用MD對鋯中輻照損傷過程進行了研究，模擬中，選定一原子賦予一定的初速度來模擬輻照過程。原子被賦予的速度與PKA能量的轉換關系為： E=mv2 (3.3)式中，E為能量（單位：J），m為原子質(zhì)量（單位：kg），v為原子速度（單位：m/s）。本模擬實驗中，原子速度v單位為À/ps，單位換算涉及常量有1ev=1.6X10-19J;C12=1.66X10-

27、27kg;1À/ps=100m/s.3.2模擬設置3.2.1 模擬體系范圍 在模擬過程中,雖然缺陷不直接與熱浴區(qū)域發(fā)生相互作用,但模擬體系中間活動區(qū)域的大小對缺陷的產(chǎn)生也有影響,這主要是因為系統(tǒng)的能量密度是由PKA能量的大小以及活動區(qū)域的體積來控制的13。3.2.2 平衡時間MD模擬的優(yōu)點在于計算量小、計算時間短。模擬過程的計算時間長短由平衡時間來決定。平衡時間是分子動力學模擬中一個重要的參數(shù),如果過大則會增加不必要的計算量；如果過小,模擬體系未達到平衡就輸出模擬結果,使模擬情況不能反映真實情況。所以選擇合適的平衡時間對于計算時間節(jié)省和,實驗結果的準確性是重要的。本文對輻照模擬的實驗

28、時間均>10ps。模擬體系的大小對平衡平衡時間也有影響。在MD模擬中,原子不斷的運動,則原子與原子的之間的相互作用力會不斷的改變,這就是能量產(chǎn)生波動的原因。隨著體系的變小,影響原子運動的條件增多,所以原子達到平衡態(tài)所需要的時間也就越長。3.2.3 弛豫時間 模擬中，需要初始化每個原子的初始速度和坐標，但初始化的矢量不一定就是對應勢函數(shù)的最小值，需要最小化來弛豫應力。體系的初始坐標和初始速度設置以后,在進行模擬體系的性質(zhì)以前,首先必須使體系進行趨于平衡的過程。在這個過程中,動能、勢能相互轉化,當動能、勢能、總能量只在某平均值附近波動、趨于平穩(wěn)時,就可認為體系就達到了平衡。3.3模擬過程和結

29、果本文以HCP密排六方晶體結構的金屬鋯（Zr）為研究對象，使用分子動力學方法對其級聯(lián)碰撞的過程進行了模擬研究。首先對鋯的結構和性質(zhì)進行簡要介紹，再結合輻照損傷理論和MD模擬方法，對鋯級聯(lián)過程的PKA入射能量的影響進行了模擬實驗，并對結果進行了分析。3.3.1 Zr的基本性質(zhì)在1789年，德國化學家Martin Klaproth分析了一塊鋯石，并分離出了鋯的“泥土”形式氧化鋯，也就是氧化物ZrO2。隨著錯合金在原子能工業(yè)上的應用,錯工業(yè)有了迅速的發(fā)展，成為一種重要的戰(zhàn)略材料,被譽為“原子時代的第一金屬”14。Zr在化學元素周期表中位于第5周期,其位置在IVB族的欽(Ti)之下,原子序數(shù)為40,平

30、均相對原子質(zhì)量為91.22。純鋯在室溫下具有HCP密排六方的晶體結構,稱為-Zr晶格常數(shù)為a=0.323nm,c=0.515nm。a-Zr在865時發(fā)生同素異晶轉變,轉變?yōu)轶w心立方-Zr,相變時會有較大的體積收縮。鋯是鈦的同族元素，也是同素異構體。862以下呈密排六方晶格，稱為鋯。862以上呈體心立方晶格，稱為鋯。純鋯的熔點Tm是1850左右。a-Zr的結構如圖3-2所示。圖3-2 -Zr的晶格結構由于核性能優(yōu)異、力學性能適中和加工性能良好,因而錯已被普遍用作核動力水冷反應堆的燃料包殼管和結構材料,如壓力管、容器管、孔道館、導向管、定位格架、端塞和其他結構文件,這是錯材的主要用途,占整個錯加工

31、材的80%15在改善鋯合金機械性能和抗腐蝕性能，研究人員已經(jīng)開發(fā)出了N18,N36等新鋯合金，其焊接性能和力學性能也有所提高。3.3.2 級聯(lián)碰撞的演化過程 100K溫度下金屬Zr中能量為1KeV的PKA原子入射模擬的結果表明：Zr的級聯(lián)碰撞過程非常短暫，大約是ps量級（1ps=10-12s）。隨著時間的推移，一開始缺陷迅速增多，隨后漸漸下降，趨于平緩，最后達到一個穩(wěn)定狀態(tài)。這是因為，當高于Zr原子離位閾值Ed的粒子入射時，晶體內(nèi)會產(chǎn)生級聯(lián)碰撞過程，產(chǎn)生間隙原子和空位并存的缺陷。在碰撞末期，所有原子的能量都降到離位閾值以下，不再撞擊晶格原子。之后，通過間隙原子的彈性復合，出現(xiàn)“急增-衰減-平緩

32、”的現(xiàn)象。利用lammps軟件對輻照損傷過程進行計算機計算，以便更直接地得到級聯(lián)碰撞過程的演化情況。觀察的結果也和缺陷產(chǎn)生規(guī)律一致。當PKA原子能量增大時，即輻照增強時，級聯(lián)碰撞將更為劇烈。這樣，將有更多的原子被賦予大于離位閾值的能量，從而有更多的原子從點陣移出。本文模擬了PKA能量分別是1KeV、5 KeV、10 KeV入射Zr基體。結果發(fā)現(xiàn)，隨著PKA能量的增大，Zr基體內(nèi)產(chǎn)生的缺陷數(shù)目越多，這說明PKA能量的增加會使碰撞加劇，即輻照損傷加劇。3.3.3 MD元胞設置如下圖所示，a-Zr為HCP結構，每個晶胞中有兩個原子，位于（0，0，0）和（1/3，2/3，1/2）處。把HCP晶格變成如

33、圖中紅線框所示的長方體晶格, 每個晶格中有4個原子，位于（0,0,0）、（0.5,5/6,0.5）、（0,1/3,0.5）和（0.5,0.5,0）。設原始坐標下晶格參數(shù)為a0和c0，則本文所采用的長方體晶格尺寸為： a0=0.3231nm, c0=5.148nm。a= a0,b=a0c=c0圖3-3 HCP結構平面圖3.3.4 模擬步驟（1）在溫度設定為600K且不變的情況下，設定一個大小為203020總原子個數(shù)為16000個原子的的長方體元胞。初始晶胞大小為a0=0.3231nm, c0=5.148nm。在600K和零外壓下運行30000個時間步長后，得到平衡狀態(tài)下的元胞體積和晶胞的參數(shù)，結

34、果為 a = 0.3232 nm b = 0.5598 nm c = 0.5212 nm(2)在不同的PKA能量下，選取不同大小的MD盒子，其選取的MD盒子表3-1:表3-1 MD盒子設計PKA能量 MD元胞大小 總原子個數(shù) 1KeV 301520 36000 5KeV 502530 150000 10KeV 603040 288000先在NV條件下運行30000步，其時間步長為1fs.。運行結束經(jīng)過馳豫過程后得到無缺陷狀態(tài)下的平衡構型（3）選取MD胞中的一個原子，賦予它不同的初始動能（由 E=mv2 可知給予原子一定的能量它就擁有一定的初始速度）,速度的方向隨機(可采用01中的隨機數(shù)得到設定

35、兩個隨機數(shù)sin，sin ,從而確定速度的方向)。在NVE下運行，采用變時間步長法，在不同能量下重復多次，最終得到平均后的實驗結果。能量和次數(shù)如表3-2所示:表3-2 不同PKA能量下的重復次數(shù) PKA能量 重復次數(shù) 1KeV 5 5KeV 10 10KeV 153.4結果分析與討論本文以密排六方晶體結構的Zr為研究對象，使用分子動力學計算模擬方法，對輻照級聯(lián)進行了模擬研究。通過計算，對模擬實施過程中的時間步長、熱裕、缺陷判別等方面進行了測試，建立了一套模擬參數(shù)。利用分子動力學計算軟件lammps對級聯(lián)過程進行了觀察。3.4.1 模擬過程中的溫升效應 NVE條件下，系綜的粒子數(shù)和總能量保持恒定

36、，系綜與外界沒有能量交換。模擬設定的溫度是600K。當PKA原子獲得模擬給予的初始能量，原子即以高速度發(fā)射，獲得動能。撞擊過程必然導致其他晶格原子的移位，造成混亂，同時由于撞擊PKA原子動能衰減。因為能量的轉化，伴隨著PKA原子動能的下降，原子的勢能會上升，但動能并非全部轉化為勢能，一部分機械能的損失轉化成了熱量。微觀世界里的劇烈撞擊，短時間內(nèi)轉化了大量熱能，這樣體系溫度迅速攀升。體系的平均溫度可達到800K。由于PKA原子的撞擊，離位原子產(chǎn)生了速度，離位原子的動能提高，再者，一部分動能的貢獻來自于高溫，即800K體系熱量的傳送。這樣，體系溫度又開始下降，最后，溫度逐漸平穩(wěn)下來。所以，整個過程

37、呈現(xiàn)一個溫度快速上升繼而逐漸下降最后趨于穩(wěn)定狀態(tài)。碰撞后體系中溫度隨時間的變化關系如圖3-43-6。圖3-4（1Kev）溫度隨時間的變化圖3-5（5Kev）溫度隨時間的變化圖3-6（10Kev）溫度隨時間的變化原子碰撞后動能轉換為勢能，動能減少，溫度下降。后因勢能弛豫，溫度小幅上升，最后處于穩(wěn)定狀態(tài)。與此同時，動能和勢能的變化如圖3-73-9：圖3-7 1Kev下動能和勢能對時間的變化關系圖圖3-8 5Kev時動能和勢能與時間的變化關系圖3-9 10Kev下動能和勢能與時間的變化關系3.4.2 級聯(lián)過程中空位與離位原子在不同時刻的分布模型圖3-10直觀展示了整個級聯(lián)過程中空位和離位原子分布的變

38、化。圖中，紅色圓點表示空位，綠色圓點表示離位原子，如圖3-10（a）（e），在PKA能量為1Kev條件下，級聯(lián)開始時，空位和離位原子都非常少，然后空位和離位原子都隨著時間增多，大約進行到0.38ps時，空位和離位原子的數(shù)量都達到峰值，如圖3-10（c）所示。之后空位和離位原子數(shù)量開始減少，最后達到一個穩(wěn)定狀態(tài)，如圖3-10（e）所示。3-11和3-12的各圖表現(xiàn)了相似的過程，不同之處在于PKA能量為5Kev和10Kev時產(chǎn)生空位和離位原子的峰值更大。 （a） （b） （c） （d） （e）圖3-10 PKA能量為1Kev時各時刻的空位與離位原子分布(a) time=ps (b) time=0.

39、18ps (c) time=0.38ps (d) time=1.37ps (e) time=3.77ps （a） （b） （c） （d） (e)圖3-11 PKA能量為5Kev時各時刻的空位與離位原子分布(a) time=ps (b) time=0.83ps (c) time=2.06ps (d) time=4.54ps (e) time=8.21ps （a） （b） （c） （d） (e)圖3-12 PKA能量為10Kev時各時刻的空位與離位原子分布(a) time=ps (b) time=0.15ps (c) time=0.76ps (d) time=3.60ps (e) time=11.

40、01ps3.4.3 級聯(lián)過程中產(chǎn)生的空位與離位原子統(tǒng)計基于對鋯的基本性質(zhì)的了解，用第2章介紹的MD模擬方法對鋯級聯(lián)的碰撞過程及此過程中設定不同PKA入射原子能量進行了模擬實驗，并對結果進行了分析。模擬結果表明：材料在輻照過程中，隨著時間的推移，晶體產(chǎn)生的缺陷數(shù)目一開始猛增，這是因為當PKA能量給出時，高于原子離位閾值的原子數(shù)目較多，碰撞概率很高。隨后，因為間隙原子與空位的彈性復合，缺陷數(shù)又會慢慢下降，最后趨于平緩，達到一個穩(wěn)定狀態(tài)；隨著PKA能量的增大，材料基體內(nèi)級聯(lián)碰撞產(chǎn)生的缺陷也會相應增多。本文的模擬結果較深入地了解了級聯(lián)過程及其影響因素、改造材料結構對其抗腐蝕性，機械性能等有著重要作用。

41、模擬的結果以時間為橫坐標，F(xiàn)lenkel 缺陷對數(shù)目為縱坐標，繪制了不同PKA入射原子能量下的關系圖。如圖3-133-15。 圖3-13 time-Flenkel 缺陷對（1Kev） 圖3-14 time-Flenkel 缺陷對（5Kev） 圖3-15 time-Flenkel 缺陷對（10Kev）3.4.4 空位和離位原子的配位數(shù)統(tǒng)計級聯(lián)發(fā)生開始前，體系中空位和離位原子很少，且都小于2，基本上都是離散的。說明發(fā)生級聯(lián)碰撞前體系中的缺陷很少，而且沒有原子團簇。大約在ps時，空位的配位數(shù)開始增大，同時有1配位的空位出現(xiàn)；大約在ps時，出現(xiàn)了2配位的空位；大約在ps時，出現(xiàn)了配位數(shù)大于2的空位。說

42、明此時已出現(xiàn)大量原子團簇。離位原子配位數(shù)的分布情況與空位相似，在此不再贅述。圖3-16ae為在10Kev條件下空位配位數(shù)不同時刻的分布情況。（a）t=0.028ps, (b)t=0.061ps, (c)t=0.082ps, (d) t=0.101ps, (e)t=0.143ps(a) (b)(c)(d)(e) 圖3-16 PKA能量為10Kev時各時刻的空位的分布4 總結本文利用分子動力學計算模擬方法，借助lammps計算機計算軟件，以HCP密排六方Zr為研究對象，對材料輻照損傷過程進行了模擬，在原子尺度上研究了目前無法用實驗的方法觀察到的輻照損傷過程。重點研究了PKA入射能量對材料輻照損傷的

43、影響。此外，在本項目的研究過程中，設計了可變時間步長對溫度、Flenkel缺陷對數(shù)目、空位和離位原子等變化關系的程序，一定程度上改進和發(fā)展了此程序。主要結論如下：1）PKA入射能量會對級聯(lián)過程產(chǎn)生一定影響。在一定溫度下，分別模擬了具有不同能量的PKA（1KeV、5 KeV、10 KeV）入射時的輻照過程。發(fā)現(xiàn)PKA能量增大時，輻照級聯(lián)過程將更劇烈。2）材料輻照作用過程非常短暫，大約是ps量級；隨時間的推進，缺陷數(shù)目迅速到達至高點，然后下降，最后趨于平緩；3）級聯(lián)發(fā)生的一開始，由于碰撞原子對其他原子進行了碰撞，溫度迅速升高。之后，動能急劇下降，離位原子與空位的逐漸產(chǎn)生使勢能得到積累，同時溫度下降，然后由于勢能弛豫以及體系內(nèi)部的彈性回復，部分離位原子與空位結合，勢能減少，溫度小幅上升，之后溫度穩(wěn)定在一定范圍內(nèi)波動。最后，溫度穩(wěn)定在720左右。過程中，溫度呈現(xiàn)一種突升急降小幅上升平緩的態(tài)勢。4）在PKA入射時體系內(nèi)缺陷很少，隨著時間的推移，