分子動力學(xué)方法模擬Be高溫和高壓的狀態(tài)方程和力學(xué)性質(zhì)

黃曉玉, 何朝政, 宋海珍, 殷復(fù)榮, 于榮梅

(南陽師范學(xué)院物理與電子工程學(xué)院, 河南 南陽 473061)

分子動力學(xué)方法模擬Be高溫和高壓的狀態(tài)方程和力學(xué)性質(zhì)

黃曉玉, 何朝政, 宋海珍, 殷復(fù)榮, 于榮梅

(南陽師范學(xué)院物理與電子工程學(xué)院, 河南 南陽 473061)

利用分子動力學(xué)方法, 采用修正鑲嵌原子勢函數(shù), 在兩組勢參數(shù)下模擬Be單晶的等溫、等壓狀態(tài)方程.通過徑向分布函數(shù), 均方根位移隨時間的變化分析Be的結(jié)構(gòu)信息.比較兩組參數(shù)的模擬結(jié)果, 第一組勢參數(shù)作用下體系原子間作用較小, 較易壓縮.兩組勢參數(shù)300K等溫壓率比較,第二組勢參數(shù)的模擬結(jié)果與實驗符合的更好.對于彈性常數(shù),第二組勢參數(shù)不能考慮溫度和壓強的影響.用第一組勢參數(shù)計算Be單晶的彈性常數(shù), 根據(jù)彈性常數(shù)隨壓強的變化判斷材料的低溫相變, 相變壓強為320GPa, 與第一性原理計算符合較好; 彈性常數(shù)隨溫度增加而減小, 與實驗結(jié)果相符.總體而言, 第一組勢函數(shù)能更全面的反映Be單晶的結(jié)構(gòu)特點.

分子動力學(xué)方法; 修正型嵌入原子勢; 狀態(tài)方程; 彈性常數(shù)

金屬鈹(Be)是一種重要的戰(zhàn)略材料, 在武器系統(tǒng)、原子能、航空航天等工業(yè)領(lǐng)域中有著很廣泛的應(yīng)用[1]. 雖然Be原子的電子排布(1s22s2)非常簡單, 但是由Be 原子組成的金屬卻具有非常獨特的物理性質(zhì)如Debye溫度高、熔點高、熔化潛熱大、熱中子吸收截面小、泊松比小等.Be金屬薄層被用作激光聚變靶丸的燒蝕層, 另外, 金屬Be還被用作核反應(yīng)堆的中子減速劑、反射體材料和中子源等, 這些都需要在高溫高壓的條件下進(jìn)行, 高溫和高壓下金屬Be結(jié)構(gòu)和性質(zhì)的研究吸引眾多研究者的注意.

近年來人們主要研究了金屬Be的相變以及彈性性質(zhì)[2-4].在室溫和一個大氣壓下金屬Be是六角密排結(jié)構(gòu)(hcp), 研究發(fā)現(xiàn)無論是改變金屬Be所處的壓強還是溫度條件都會使金屬Be的結(jié)構(gòu)由hcp轉(zhuǎn)變?yōu)轶w心立方結(jié)構(gòu)(bcc).常壓(保持一個大氣壓下)下, 改變溫度發(fā)現(xiàn)金屬Be在達(dá)到熔化之前大約30K時就由hcp相轉(zhuǎn)變?yōu)閎cc相[5-6].在實驗方面, 使用金剛石壓砧(DAC) 壓縮金屬Be, 通過觀察電阻的變化來測量金屬Be相變的壓強, 一直到40 GPa時, Reichlin等人[7]還是一直沒有觀察到相變, 金屬Be仍然是 hcp 結(jié)構(gòu).進(jìn)一步提高實驗手段, 采用大功率的X 射線衍射方法觀測, 還是采用DAC壓縮金屬Be直到壓強達(dá)到171GPa, Nakano 等人[8]還是沒有觀察到晶體Be結(jié)構(gòu)的變化.目前更高壓強下的實驗數(shù)據(jù)還沒有找到.在理論方面, Benedict等人[9]利用第一性原理方法, 發(fā)現(xiàn)金屬Be在400GPa附近時才發(fā)生相變, 并做出了Be的hcp相、bcc相和液相的三相圖, 但是目前實驗上還得不到驗證.Fen Luo等人[10]用第一性原理方法計算高溫、高壓下金屬Be的熱力學(xué)性質(zhì), 給出相變壓強為390Gpa, 其計算結(jié)果中高壓下的彈性常數(shù)C44較實驗數(shù)據(jù)偏大.Aimin Hao等人[11]用第一性原理方法計算高壓下金屬Be的結(jié)構(gòu)穩(wěn)定性和彈性常數(shù), 給出相變壓強為320Gpa.對與彈性常數(shù), Marie等人[12]利用超聲方法研究金屬Be的彈性模量在300K到1000K之間隨溫度的變化趨勢, 數(shù)據(jù)和原有實驗有明顯的差異.K.Kádas等人[2]利用第一性原理計算各彈性常數(shù), 體積彈性模量隨溫度的變化.總之, 關(guān)于金屬Be相變的壓強, 第一性原理相關(guān)研究結(jié)果還存在較大的爭議.另外, 金屬Be彈性模量隨溫度的變化與實驗結(jié)果差別較大.

目前, 對于金屬Be的高溫高壓下的熱力學(xué)性質(zhì)的理論研究還不是很詳盡.因此, 我們需要采用其他方法和模型來進(jìn)行研究.本文采用分子動力學(xué)的方法和和兩組修正型嵌入原子勢(Modified embedded atom method (MEAM))參數(shù)計算高溫、高壓下金屬Be的狀態(tài)方程和彈性常數(shù), 將計算數(shù)據(jù)與實驗和第一性原理的結(jié)果進(jìn)行比較, 并給出相關(guān)分析.

1 計算方法

分子動力學(xué)方法是基于牛頓力學(xué)基本原理, 通過對運動方程積分得到系統(tǒng)粒子的位臵構(gòu)型和動量, 并運用統(tǒng)計力學(xué)方法得到物質(zhì)的宏觀性質(zhì).本文模擬體系是晶格結(jié)構(gòu)為hcp結(jié)構(gòu), 晶格常數(shù)a為2.27?, c為3.59?, 含有2048個Be原子的塊金屬固體.初始速度是按300K溫度下麥克斯韋波爾茲曼分布提取的, 每次模擬的結(jié)果都是在系統(tǒng)達(dá)到平衡之后進(jìn)行統(tǒng)計平均的.模擬體系選取的時間步長是1fs, 選擇三維周期性邊界條件, 選用的系綜包括等溫等壓系綜(NPT)、和正則系綜(NVT)等.

分子動力學(xué)(MD)方法是利用原子間相互作用勢來計算作用在原子上的力.因此, 相互作用勢就是模擬最關(guān)鍵的輸入要素.本文中, 采用Baskes提出的hcp結(jié)構(gòu)金屬的MEAM相互作用勢[13].由Daw和Baskes提出的EAM勢, 能夠很好的描述金屬間的相互作用[14-15].他們從準(zhǔn)原子理論出發(fā), 把系統(tǒng)中的每一個原子均看作是嵌在由其它原子形成的主晶格中的雜質(zhì).運用密度泛函理論并考慮局域密度近似得到系統(tǒng)的總能量為:

我們有必要對本文用到的兩組勢參數(shù)給出說明和解釋.其中, 第二組勢參數(shù)是1993年M I Baskes等人[13]以HCP結(jié)構(gòu)金屬為研究對象開發(fā)的多體勢, 這組勢參數(shù)只有在溫度為0K的情況下才能很好的模擬材料的力學(xué)、熱力學(xué)性質(zhì), 即使在很小的溫度下Be的晶體結(jié)構(gòu)就會變得很不穩(wěn)定.第一組勢參數(shù)是2007年V.V.Dremov等人[16]在MEAM理論框架下, 對勢參數(shù)進(jìn)行新的優(yōu)化得到的.優(yōu)化后的參數(shù)能夠很好的模擬材料的力學(xué), 熱力學(xué)性質(zhì).下面我們用表格的形式分別給出兩組勢參數(shù):

表1 金屬Be hcp相的MEAM參數(shù)其中）是結(jié)合能,是平衡時的最近鄰距離,是能量函數(shù)的指數(shù)衰退因子,是嵌入函數(shù)的屏蔽因子,是電子密度的指數(shù)衰退因子,t（h）是電子密度的權(quán)重因子.Table 1 Parameters of the MEAM potential for beryllium(,rc(?)equilibrium nearest-neighbor distance,the exponential decay factor for the universal energy function,A the scaling factor for the embedding function,the exponential decay factor for the atom densities,t(h)the weighting factors )

表1 金屬Be hcp相的MEAM參數(shù)其中）是結(jié)合能,是平衡時的最近鄰距離,是能量函數(shù)的指數(shù)衰退因子,是嵌入函數(shù)的屏蔽因子,是電子密度的指數(shù)衰退因子,t（h）是電子密度的權(quán)重因子.Table 1 Parameters of the MEAM potential for beryllium(,rc(?)equilibrium nearest-neighbor distance,the exponential decay factor for the universal energy function,A the scaling factor for the embedding function,the exponential decay factor for the atom densities,t(h)the weighting factors )

2 結(jié)果和討論

結(jié)果與討論部分由三個部分組成: 1)比較兩組勢參數(shù)下的高溫、高壓狀態(tài)方程, 了解兩組勢參數(shù)的特點和區(qū)別; 2)比較兩組勢參數(shù)下的徑向分布函數(shù), 均方位移隨時間的變化, 了解兩組勢參數(shù)的特點和區(qū)別; 3.)用第一組勢參數(shù)計算體系高溫、高壓下的彈性常數(shù), 體積彈性模量, 將結(jié)果與對應(yīng)的第一性原理計算、實驗結(jié)果比較.

2.1 體系的狀態(tài)方程

首先, 利用兩組勢函數(shù)模擬體系的等溫狀態(tài)方程曲線(恒溫300K), 如圖1所示: 在第一組勢函數(shù)作用下, 平衡后體系的體積V1, 相應(yīng)的第二組勢函數(shù)作用下體系的體積V2.壓強小于150Gpa時:V1>V2,壓強增加至150Gpa時:V1=V2, 壓強大于150Gpa時, V1

圖1 兩組勢參數(shù)300K等溫狀態(tài)方程比較Fig.1 Isothermal equation of state at 300K with two potential parameters

圖2 兩組勢參數(shù)300K等溫壓率比較 Fig.2 Isothermal compressibility at 300K with two potential parameters.

2.2 體系的結(jié)構(gòu)分析

在統(tǒng)計力學(xué)中, 徑向分布函數(shù)(rdf)指的是給定某個粒子的坐標(biāo), 其他粒子在空間的分布幾率.徑向分布函數(shù)可以用來研究物質(zhì)的有序性, 能夠很好的描述非晶和晶體材料的原子結(jié)構(gòu)特征.如圖4所示: 溫度為300K、壓強為0.1Mpa時, 兩種勢函數(shù)作用下體系的徑向分布函數(shù).在第一組勢函數(shù)作用下, 三個峰的位臵分別對應(yīng)第一、二、三最近鄰距離: 2.25?,4.15?,6.15?, 第二組勢參數(shù)相應(yīng)的值分別為: 2.25?,3.75?,5.65?.可以看出第一組勢函數(shù)作用下, 體系長程有序結(jié)構(gòu)更松弛.

分子動力學(xué)計算系統(tǒng)中的原子由起始位臵不停移動, 每一瞬間各原子的位臵皆不相同.粒子位移平方的平均值稱為均方位移(MSD).如圖5a,b)所示: 體系在第一組勢函數(shù)作用下,常壓、溫度分別為300K和1500K時均方根位移隨時間的變化.300K時, MSD隨時間增加變化不大, 1500K時MSD隨時間增加線性增長; 如圖6a,b)所示: 體系在第二組勢函數(shù)作用下常壓、溫度分別為300K和1500K時均方根位移隨時間的變化.同樣, 300K時, MSD隨時間增加變化不大, 1500K時MSD隨時間增加線性增長.這說明1500K時, Be晶體的原子結(jié)構(gòu)已趨于無序, 基本熔化, 與已知的熔點1550K基本一致[18].兩圖的區(qū)別在于: 第一組勢函數(shù)作用下體系的MSD數(shù)值較大; 1500K時第一組勢函數(shù)作用下MSD隨時間增加線性增長的斜率較大, 易膨脹.

圖3 兩種勢函數(shù)等壓狀態(tài)方程比較Fig.3 Isobaric equation of state with two potential parameters

圖4 兩種勢函數(shù)徑向分布函數(shù)比較 Fig.4 Radial distribution function at 300K with two potential parameters

圖5 勢參數(shù)1下體系的均方位移Fig.5 Mean square displacement with the first potential parameters

圖6 勢參數(shù)2下體系的均方位移 Fig.6 Mean square displacement with the second potential parameters

表2 兩種勢函數(shù)下等壓狀態(tài)方程分段擬合Table 2 Curve fitting of isobaric equation of state

2.3 體系的彈性常數(shù)

材料的彈性常數(shù)描述了它對所加外力的響應(yīng), 或者可以這樣說, 彈性常數(shù)描述了為維持一個給定的形變所需的應(yīng)力.對體系施加不同的變形, 根據(jù)總能量和應(yīng)變的關(guān)系, 擬合彈性常數(shù)的具體數(shù)值.其中第二組勢參數(shù)的計算結(jié)果有大的無規(guī)則波動, 不能有效考慮溫度和壓強的效應(yīng).因此, 將第一組勢參數(shù)模擬的彈性常數(shù)和相應(yīng)的實驗數(shù)據(jù)、第一性計算數(shù)據(jù)比較, 給出分析結(jié)果.

圖7 MD模擬各向異性彈性常數(shù)值Fig.7 Anisotropic elastic constants by MD simulation

圖8 MD模擬單軸彈性模量 Fig.8 Uniaxial elastic modulus by MD simulation

圖9 體積彈性模量隨壓強變化（MD與第一性原理比較）Fig.9 The variation of Bulk modulus with pressure (MD simulation and first principles study )

圖10 MD體積彈性模量值與實驗值、第一性計算值比較 Fig.10 The variation of Bulk modulus with temperature (experiment and first principles study )

對于HCP結(jié)構(gòu)的Be晶體, 它的獨立的彈性常數(shù)為: C11, C12, C13, C33和C44.如圖7所示: 壓強范圍為0～300Gpa各彈性常數(shù)隨壓強的變化.在指定壓強范圍內(nèi), 各彈性常數(shù)隨壓強的變化線性增加.其中, C11和C33數(shù)值接近, 變化趨勢很接近, 隨壓強增加速度最快; C44隨壓強增加速度最慢; C12和C13數(shù)值接近, 變化趨勢也很接近, 隨壓強增加速度居中.不同彈性常數(shù)的增加并不是完全等比例的, 這說明, Be單晶的彈性各向異性關(guān)系在不同壓強下也不一樣.HCP結(jié)構(gòu)的Be晶體的單軸彈性常數(shù)分別為: Ca, Cb, Cc.Ca=0.5*(C11-C12), Cb=C44, Cc=1/6*(C11+2*C33+C12-4*C13), 如圖8所示: 壓強范圍為0～300Gpa, 單軸彈性常數(shù)隨壓強的變化.其中, Cb,Cc的數(shù)值隨壓強增加速度較快; Cc的數(shù)值隨壓強增加速度較慢; 壓強在150Gpa～300Gpa范圍, Ca的數(shù)值隨壓強增加變化趨于平緩, 數(shù)值變化很小.

如前所述, Be單晶的常溫高壓相變行為: hcp～bcc, 相變壓強的確定尚存爭議.本文通過Be單晶體積彈性模量隨壓強的變化給出分子動力學(xué)方法對于相變的判斷.如圖9所示, 給出壓強范圍為0～400Gpa, MD方法和第一性原理計算方法, Be單晶的體積彈性模量隨壓強的變化.在0～300Gpa范圍內(nèi), MD模擬結(jié)果與Hao A等人的計算結(jié)果符合的很好, 和Luo F的數(shù)據(jù)有較小的差別. 壓強為320Gpa時, 體積彈性模量發(fā)生突變, 可以推測Be單晶發(fā)生相變行為.這一結(jié)果和Hao A等人給出的相變壓強一致, 和Luo F等人給的結(jié)果有一定偏差.如圖10所示: 三種方法體積彈性模量隨溫度變化, 彈性常數(shù)隨溫度的增加而減小.MD計算結(jié)果與實驗數(shù)據(jù)符合的較好[2],第一性原理方法計算數(shù)據(jù)與實驗數(shù)據(jù)相比有較大偏差.

3 結(jié)論

本研究利用兩組勢函數(shù)模擬Be單晶的等溫、等壓狀態(tài)方程; 徑向分布函數(shù)和均方根位移; 彈性常數(shù)、體積彈性模量.可以得出結(jié)論: 第一組勢函數(shù)作用下, 體系的可壓縮性強, 原子間的作用力較弱; 長程有序結(jié)構(gòu)更加松弛; 均方位移數(shù)值較大且隨時間的變化率值較大.

兩組勢參數(shù)300K等溫壓率比較,第二組勢參數(shù)的模擬結(jié)果與實驗符合的更好.但是, 第二組勢函數(shù)計算彈性常數(shù)時, 不能有效地考慮壓強和溫度的效應(yīng).第一組勢函數(shù)模擬的彈性常數(shù)隨壓強的變化顯示體系的各項異性特點, 其中C軸可壓縮性較強; 體積彈性模量在壓強為320Gpa時有突變, 據(jù)此可推測Be單晶hcp到bcc結(jié)構(gòu)的相變行為, 這一結(jié)果和已有的第一性原理計算符合較好.體積彈性模量隨溫度增加而減小, 與已有的實驗數(shù)據(jù)符合較好, 彌補第一性原理數(shù)據(jù)的偏差.本文將兩組勢函數(shù)的模擬結(jié)果比較分析, 研究結(jié)果說明: 第一組勢函數(shù)能夠更全面的反映Be單晶的結(jié)構(gòu)特點.準(zhǔn)確的勢函數(shù)是分子動力學(xué)模擬的關(guān)鍵要素, 希望我們的研究結(jié)果對后來的理論模擬工作提供參考.

[1]張友壽, 秦有鈞, 吳東周, 等.鈹?shù)姆勰┮苯鸸に嚰昂附友芯窟M(jìn)展[J].焊接學(xué)報, 2001, 22: 6-10.

[2]KáDAS K, VITOS L, AHUJA R, JOHANSSON B, et al.Temperature-dependent elastic properties of α-beryllium from first principles[J].Phys Rev B, 2007, 76: 235109.

[3]ROBERT G, LEGRAND P, BERNARD S.Multiphase equation of state and elastic moduli of solid beryllium from first principles[J].Phys Rev B, 2010, 82:104118.

[4]宋海峰, 劉海風(fēng).金屬鈹熱力學(xué)性質(zhì)的理論研究[J].物理學(xué)報, 2007, 56:2833-2838.

[5]MARTIN A J, MOORE A.1959 J.Less-Common Met.1 85.

[6]FRANCOIS M M.contre, Proc Int Conf Grenoble.Presses Univ de France, Paris, 1965.

[7]REICHLIN R L.Measuring the electrical resistance of metals to 40 GPa in the diamond-anvil cell[J].Rev Sci Instrum, 1983, 54: 1647-1652

[8]NAKANO K, AKAHAMA Y, KAWAMURA H.X-ray diffraction study of Be to megabar pressure[J].J Phys: Condens Matter, 2002, 14: 10569-10575.

[9]BENEDICT L X, OGITEU T, TRAVE A, et al.Calculations of high pressure properties of beryllium: Construction of a multiphase equation of state[J].Phys Rev B, 2009, 79: 064106-064111.

[10]LUO F, CAI L C, CHEN X R, JING F Q, et al.Ab initio calculation of lattice dynamics and thermodynamic properties of beryllium[J].J Appl Phys, 2012, 111: 053503-053508.

[11]HAO A, ZHU Y.First-principle investigations of structural stability of beryllium under high pressure [J].J Appl Phys, 2012, 112: 023519-023524.

[12]NADAL M H,BOURGEOIS L.Elastic moduli of beryllium versus temperature: Experimental data updating [J].J Appl Phys, 2010, 108: 033512-033517.

[13]BASKES M I, JOHNSON R A.Modified embedded atom potentials for HCP metals[J].Modelling Simul Mater Sci Eng, 1994, 2: 147-152.

[14]DAW M S, BASKES M I.Embedded-atom method: Derivation and application to impurities, surfaces, and other defects in metals[J].Phys Rev B, 1984, 29: 6443-6448.

[15]DAW M S, FOILES S M, BASKES M I.The embedded-atom method: a review of theory and applications[J].Mater Sci Rep, 1993, 9: 251-256.

[16]DREMOV V V, KARAVAEV A V, Kutepov A L, et al.Molecular Dynamics Simulation of Thermodynamic and Mechanical Properties of Be and Mg[J].AIP Conf Proc, 2007, 955: 305-310.

[17]VELISAVLIEVIC N, CHESNUT G N, VOHRA Y K.Structural and electrical properties of beryllium metal to 66 GPa studied using designer diamond anvils[J].Phys Rev B, 2002, 65: 172107-172112.

Study of equations of state and mechanical property of Be by molecular dynamics method

HUANG Xiao-yu, HE Chao-zheng, SONG Hai-zhen, YIN Fu-rong, YU Rong-mei
( School of Physics and Electronic Engineering, Nanyang Normal University, Nanyang 473061, P.R.C.)

Beryllium is an important strategic nuclear material, which was used in weapons system, atomic energy, aerospace fields. Based on these applying requirements, study of beryllium on structure and properties was necessary. This paper uses molecular dynamics method and modified embedded atom potential to study equations of state of beryllium under high temperature and high pressure with two sets of potential parameters. Structure information of beryllium was analyzed by the radial distribution function and mean square displacement. To compare the computer simulation results with two sets of potential parameters, with the first potential parameters, forces between atoms were smaller, more compressible. For the elastic constants, the second potential parameters can not consider the effects of temperature and pressure; elastic constants changed with temperature and pressure were computed with the first potential parameters. The results were in good accordance with the experiments and the first principles calculations. There is a sudden change when the elastic constant changed with the pressure which can be guessed as the phase transition. The pressure was 320Gpa which was in accordance with the the first principles calculations. In short, the first potential parameters can reflect the structural characteristics of beryllium more comprehensive than the first one.

molecular dynamics method; modified embedded atom potential; equation of state; elastic constant

O52

A

1003-4271(2014)06-0921-07

10.3969/j.issn.1003-4271.2014.06.21

2014-10-08

黃曉玉(1985-), 女, 漢族, 河南南陽人, 博士, 講師, 研究方向: 材料狀態(tài)方程、熱力學(xué)性質(zhì), E-mail: huangxiaoyu198510@163.com.

國家自然科學(xué)基金青年科學(xué)基金(批準(zhǔn)號: 11304167); 河南省教師教育改革重點項目(批準(zhǔn)號: 2013-JSJYZD-053); 南陽師范學(xué)院科研基金(批準(zhǔn)號: ZX2013020、ZX2014088).