周幼華，李雯慧，陶蕾，鄭廣

（江漢大學(xué)物理與信息工程學(xué)院，湖北武漢430056）

鐵熔體快速冷卻過程的分子動(dòng)力學(xué)模擬

周幼華，李雯慧，陶蕾，鄭廣

（江漢大學(xué)物理與信息工程學(xué)院，湖北武漢430056）

利用LAMMPS建立8×8×8原子箱，采取EAM模擬大量鐵原子體系。先對(duì)體系升溫，讓體系充分均勻，然后對(duì)體系進(jìn)行不同冷卻速度下的降溫。經(jīng)過分子動(dòng)力學(xué)模擬得到能量溫度曲線、徑向分布函數(shù)（RDF）和體系結(jié)構(gòu)。結(jié)果表明：冷卻速度低于1.8×1010K/s時(shí)，結(jié)晶相變點(diǎn)和熔點(diǎn)基本重合；冷卻速度超過1012K/s時(shí)，生成物為非晶；鐵熔體的過冷度可達(dá)到700 K；冷卻速度在1011K/s附近時(shí)，能夠生成體心立方的Fe單晶。

分子動(dòng)力學(xué)模擬；快速冷卻；鐵熔體

0 引言

對(duì)于相變溫度低于熔點(diǎn)的晶體不能采用傳統(tǒng)熔融提拉方法生長(zhǎng)，因而合成單晶存在較大困難。如：傳統(tǒng)粉末冶金法在合成的FeSi2合金時(shí)趨向于形成高溫穩(wěn)定的α相，很難直接得到低溫穩(wěn)定的β相。而采用脈沖激光沉積法，由于激光引起的微液滴是主要的傳輸方式，經(jīng)歷了快速冷卻過程，直接得到β-FeSi2，而無需經(jīng)過高溫穩(wěn)定的α相［1］。與此類似，單質(zhì)鐵在不同溫度具有不同的晶體結(jié)構(gòu)，室溫下為體心立方點(diǎn)陣（α-Fe），晶格常數(shù)為0.286 3 nm；在912℃，α-Fe轉(zhuǎn)變?yōu)榫哂忻嫘牧⒎近c(diǎn)陣的γ-Fe，晶格常數(shù)為0.359 1 nm；溫度繼續(xù)升高，超過1 394℃時(shí)，將再次轉(zhuǎn)變?yōu)轶w心立方點(diǎn)陣δ-Fe；δ-Fe的體心立方結(jié)構(gòu)將一直保持到熔點(diǎn)1 538℃；在熔點(diǎn)以上，鐵由固態(tài)轉(zhuǎn)化為液態(tài)。正因?yàn)槿绱?，傳統(tǒng)的熔融提拉法很難獲得鐵單晶，通常的固態(tài)單質(zhì)鐵都是α-Fe和γ-Fe（或者是滲碳后的奧氏體）的混合體［2］。熔體中原子沒有確定的位置，結(jié)構(gòu)具有不穩(wěn)定性和不確定性，熔體的團(tuán)簇結(jié)構(gòu)，在高速冷卻時(shí)沒有足夠的時(shí)間形成晶體，成為非晶凝固液體系統(tǒng)；液相中的團(tuán)簇結(jié)構(gòu)，包括最近鄰原子的短程有序，被保持在了非晶狀態(tài)下?？焖倮鋮s過程可能會(huì)出現(xiàn)直接得到低溫穩(wěn)定相而不經(jīng)過高溫中間相。很多晶體的相變溫度低于熔點(diǎn)，因而無法采用傳統(tǒng)熔融提拉方法生長(zhǎng)，熔體快速冷卻法為這類晶體提供了一條可能的合成途徑。因而對(duì)單質(zhì)鐵熔體快速冷卻過程結(jié)晶現(xiàn)象的研究，有助于揭示這類單晶的生長(zhǎng)相規(guī)律?？焖倌讨胁⒉粷M足穩(wěn)態(tài)形核理論。實(shí)驗(yàn)利用快冷和深冷可使凝固的合金得到優(yōu)異的性能，如擴(kuò)大了固溶的極限，能使合金得到超細(xì)的晶粒度，減少偏析或者無偏析，使合金形成亞穩(wěn)相。還可能生成高的點(diǎn)缺陷密度晶體材料，并且還可以通過快速冷卻開發(fā)出準(zhǔn)晶體、非晶體、納米晶體等材料［3-5］。

本文采用分子動(dòng)力學(xué)方法模擬鐵熔體的快速冷卻過程，研究冷卻速度對(duì)熔體相變溫度的影響。通過控制模擬過程中的起始溫度（Tstart）、最終溫度（Tstop）、步長(zhǎng)（Timestep）和運(yùn)行步數(shù)（Run）來控制降溫的速度，得到每一步體系中分子的運(yùn)動(dòng)情況，計(jì)算體系能量隨溫度的變化，體系中分子的位置、動(dòng)量、能量和最近鄰原子數(shù)，由此得出鐵熔體在不同冷卻速度下的冷卻過程中呈現(xiàn)的狀態(tài)。

1 系統(tǒng)參數(shù)的計(jì)算方法

1.1 溫度

在多體體系中，溫度的定義采用在系統(tǒng)的每一個(gè)自由度上能量均分，在實(shí)際過程中，由于體系的總動(dòng)能漲落，體系的瞬時(shí)溫度也隨之改變：

其中Nf為體系自由度，m為分子質(zhì)量，v為分子速度，kB為玻爾茲曼常數(shù)，N為體系的分子數(shù)。

1.2 徑向分布函數(shù)（RDF）

RDF即原子徑向上的原子密度與體系總密度的比值函數(shù)。

其中n(r)為半徑在r到r+Δr的球殼內(nèi)的原子數(shù)，p0為理想晶體的原子密度。

1.3 嵌入原子勢(shì)（EAM）

EAM模型在金屬系統(tǒng)的研究方面取得巨大成功，因?yàn)樗軌蚝芎玫孛枋鼋饘俳Y(jié)合能的形成，類似于將金屬的原子核嵌入自由電子氣的性質(zhì)。基本思想是把晶體的總勢(shì)能分成兩部分，一部分是在晶格點(diǎn)陣上的原子核之間的相互作用，另一部分是原子核在電子云中的嵌入勢(shì)。

總勢(shì)能為

其中Ei(pi)為嵌入勢(shì)能，第二項(xiàng)是對(duì)勢(shì)項(xiàng)，pi是除第i個(gè)原子以外的所有其他原子的核外電子在第i個(gè)原子處產(chǎn)生的電子云密度之和，

2 模擬方法

采用一個(gè)8×8×8的原子箱，用EAM來計(jì)算。分子動(dòng)力學(xué)程序構(gòu)成方式：

1）讀入指定運(yùn)算條件的參數(shù)，如初始溫度、粒子數(shù)、密度、時(shí)間等。

2）體系初始化，即選定初始位置和速度。

3）計(jì)算作用于所有粒子上的力。

4）求解方程。

5）中心循環(huán)完成，計(jì)算并輸出測(cè)量值。

對(duì)于過程的模擬分兩步。先對(duì)體系升溫，讓體系充分均勻，然后對(duì)體系進(jìn)行不同冷卻速度下的降溫。通過確定Timestep、Tstart和Tstop，可以通過控制Run來控制降溫速度。

命令為

fix ID group-ID style_name keyword value...

對(duì)于style_name，采取等溫等壓系綜（NPT）。對(duì)于keyword，包括溫度標(biāo)準(zhǔn)（Tstart、Tstop、Tdamp）和壓強(qiáng)標(biāo)準(zhǔn)（Pstart、Pstop、Pdamp）。對(duì)于溫度標(biāo)準(zhǔn)，可確定Tstart和Tstop。對(duì)于壓強(qiáng)標(biāo)準(zhǔn)，可確定Pstart和Pstop。為了使計(jì)算結(jié)果能夠相對(duì)準(zhǔn)確并且防止原子丟失，關(guān)于溫度和壓強(qiáng)的弛豫時(shí)間Tdamp和Pdamp可以采取Timestep的10～100倍。

有關(guān)Timestep（單位為ps）的選取，和原子運(yùn)動(dòng)時(shí)間尺度有關(guān)。一般來說，簡(jiǎn)單原子系統(tǒng)取10 fs；剛性分子取5 fs；軟性分子，限制鍵長(zhǎng)取2 fs；軟性分子，無鍵長(zhǎng)限制取1 fs或0.5 fs。

可以通過控制Run來確定降溫速度。

采取等壓等溫系綜（NPT），步長(zhǎng)0.001ps，運(yùn)行10 000步從1 000 K升溫至2 000 K，再從2 000 K降溫至200 K，運(yùn)行1 000 000步，降溫速度為1.8×1012K/s。

3 模擬結(jié)果與分析

3.1 鐵熔體冷卻過程中的相變

模擬中得到冷卻速度從1.8×1012K/s到1.8× 1010K/s的冷卻曲線。圖1是冷卻速度為1.8×1012K/s時(shí)，體系總能量和體積與溫度的關(guān)系曲線。

圖1所示體系總能量和體積隨溫度變化曲線是連續(xù)的，基本為一條直線，沒有觀察到明顯的一級(jí)相變。因而，在此冷卻速度下鐵熔體會(huì)形成過冷液體——非晶態(tài)。

圖1 冷卻速度為1.8×1012K/s時(shí)，體積和總能量與溫度的關(guān)系Fig.1Curve of total energy and volume vs temperature at cooling rate of 1.8×1012K/s

圖2 冷卻速度為6.5×1011K/s時(shí)，體積和總能量與溫度的關(guān)系Fig.2Curve of total energy and volume vs temperature at cooling rate of 6.5×1011K/s

當(dāng)冷卻速度為6.5×1011K/s時(shí)，體系總能量和體積與溫度的關(guān)系曲線如圖2所示。能量—溫度曲線總趨勢(shì)為線性關(guān)系。出現(xiàn)了斷開，意味著相變的出現(xiàn)。當(dāng)冷卻速度為6.5×1011K/s時(shí)，相變發(fā)生在800 K；當(dāng)冷卻速度為6.0×1011K/s時(shí)，相變發(fā)生在935 K；當(dāng)冷卻速度為4.5×1011K/s時(shí)，相變發(fā)生在957 K；當(dāng)冷卻速度為3.6×1011K/s時(shí)，相變發(fā)生在1 079 K；當(dāng)冷卻速度為1.8×1011K/s時(shí)，相變發(fā)生在1 057 K；當(dāng)冷卻速度為1.2× 1011K/s時(shí)，相變發(fā)生在1 109 K；當(dāng)冷卻速度為9×1010K/s時(shí)，相變發(fā)生在1 229 K；當(dāng)冷卻速度為7.2×1010K/s時(shí)，相變發(fā)生在1 405 K；當(dāng)冷卻速度為6×1010K/s時(shí)，相變發(fā)生在1 519 K；當(dāng)冷卻速度為3.6×1010K/s時(shí)，相變發(fā)生在1 718 K；當(dāng)冷卻速度為1.8×1010K/s時(shí)，相變發(fā)生在1 810 K。

計(jì)算表明：①當(dāng)冷卻速度高于1.2×1011K/s時(shí)，固液相變發(fā)生的溫度低于912℃（1 185 K），位于體心立方的α-Fe是穩(wěn)定相；②當(dāng)冷卻速度高于6×1010K/s低于1.2×1011K/s時(shí)，固液相變溫度位于面心立方γ-Fe的穩(wěn)定區(qū)；③當(dāng)冷卻速度低于1.8×1010K/s時(shí)，固液相變溫度接近1 538℃（1 801 K），將形成過冷液體或體心立方結(jié)構(gòu)的δ-Fe。

得到的結(jié)晶溫度與冷卻速度的關(guān)系見圖3。由以上結(jié)果可知，隨著冷卻速度的降低，發(fā)生相變的溫度增大，過冷度降低。冷卻速度低于1.8× 1010K/s時(shí)結(jié)晶相變點(diǎn)和熔點(diǎn)基本重合，理論計(jì)算快速冷卻過程中，鐵熔體的過冷度可達(dá)到800K。由于冷卻速度大于1.8×1012K/s觀察不到相變，因此能形成晶體的區(qū)域過冷度可達(dá)到700 K。

3.2 鐵熔體冷卻過程中的徑向分布函數(shù)（RDF）

RDF表現(xiàn)了體系中原子分布和徑向位置的關(guān)系，對(duì)于判斷體系結(jié)構(gòu)有著重要意義。在溫度較低時(shí)，RDF的峰清晰尖銳，各個(gè)峰所對(duì)應(yīng)的值，對(duì)應(yīng)原子周圍最近鄰、次近鄰配位的位置，體現(xiàn)了規(guī)律的周期性晶體結(jié)構(gòu)；隨著溫度升高，晶體的周期性結(jié)構(gòu)逐漸解體，RDF的峰變寬，一些峰消失；直至到達(dá)液相，RDF的值不表示配位情況，而是反映了此時(shí)體系原子相對(duì)于中心原子位置的概率。

體系以6.5×1011K/s的速度降溫時(shí)，末態(tài)結(jié)構(gòu)的RDF如圖4所示，結(jié)晶后2近鄰原子、3近鄰原子和4近鄰原子最可幾分布點(diǎn)，分別在1.15、1.65和1.98倍的最近鄰位置，這一結(jié)論和體心立方的格點(diǎn)排布相同。

圖3 結(jié)晶溫度與冷卻速度的關(guān)系曲線Fig.3Curve of crystallization temperature and cooling rate

圖4 冷卻速度為6.5×1011K/s時(shí)，末態(tài)的RDFFig.4Curve of radial distribution function of iron in the final structure at cooling rate of 6.5×1011K/s

由Jackson界面平衡結(jié)構(gòu)理論模型，晶體生長(zhǎng)與相變時(shí)的溫度（TCR）和前后的熵變（ΔS）有關(guān)。Jackson常數(shù)

當(dāng)α＜2時(shí)無晶面生長(zhǎng)，2＜α＜10時(shí)在生長(zhǎng)過程中產(chǎn)生晶面。平衡情況下Fe液固相變?chǔ)?1.02［6］，因而很難獲得Fe單晶。考慮到快速冷卻過程凝固溫度TCR大幅度降低，同時(shí)相變過程中熵變?chǔ)也增大，快速冷卻過程大幅度地提高了Jackson常數(shù)α。上述結(jié)果表明，在冷卻速度為1011K/s附近可以生成體心立方Fe單晶。

4 結(jié)論

在分子動(dòng)力學(xué)模擬鐵熔體的快速冷卻過程中，通過對(duì)8×8×8個(gè)原子的體系進(jìn)行不同速度的降溫，得到能量曲線和RDF，得出如下結(jié)論：

1）冷卻速度低于1.8×1010K/s時(shí)，結(jié)晶相變點(diǎn)和熔點(diǎn)基本重合，鐵熔體的過冷度可達(dá)到700 K。

2）鐵熔體冷卻速度為1012K/s以上時(shí)，生成物為非晶。

3）在冷卻速度為1011K/s附近可以生成體心立方Fe單晶。

（References）

［1］ZHOU Y H，NIE C，TIAN H Y，et al.Investigation on the synthesis mechanism ofβ-FeSi2prepared by pulsed laser deposition［J］.Wuhan University Journal of Natural Sciences，2012，17（1）：61-66.

［2］陳景榕，李承基.金屬與合金中的固態(tài)相變［M］.北京：冶金工業(yè)出版社，1997：13-16.

［3］邵琛瑋，王振華，李艷男，等.AuCu249合金團(tuán)簇?zé)岱€(wěn)定性的原子尺度計(jì)算研究［J］.物理學(xué)報(bào)，2011，60（8）：083602-083610.

［4］王麗，邊秀房，李輝.金屬Cu熔化結(jié)晶過程的分子動(dòng)力學(xué)模擬［J］.化學(xué)物理學(xué)報(bào)，2000，13（5）：544-550.

［5］陳光，傅恒志.非平衡凝固新型金屬材料［M］.北京：科學(xué)出版社，2004：1-19.

［6］MOULSON A J，HERBERT J M.Electroceramics：ma?terials，properties，applications［M］.2nd ed.Hobo?ken，NJ：John Wiley&Sons，2003：122-125.

（責(zé)任編輯：曾婷）

Molecular Dynamics Simulation of Rapid Cooling Process of Iron Melt

ZHOU Youhua，LI Wenhui，TAO Lei，ZHENG Guang
（School of Physics and Information Engineering，Jianghan University，Wuhan 430056，Hubei，China）

Establishes 8×8×8 atomic box with LAMMPS，adopts EAM（embedded atom model）potential，to simulate large number iron atoms system.At first，heats up the system，then cool down the system under different cooling rate.Through molecular dynamics simulation，obtains energy-tem?perature curve，radial distribution function（RDF）and architecture.The results show that the crystal?lization temperature point is basically in coincident with the melting point when the cooling rate is lower than 1.8×1010K/s；the product is amorphous when the cooling rate is higher than 1012K/s；the degree of supercooling can reach 700 K；Fe single crystal can be prepared at cooling rate around 1011K/s.

molecular dynamics simulation；rapid cooling；iron melt

O561.1

A

1673-0143（2014）03-0037-04

2014-04-02

國家自然科學(xué)基金資助項(xiàng)目（61240056）

周幼華（1969—），男，副教授，博士，研究方向：功能薄膜、半導(dǎo)體薄膜、新型光電器件的開發(fā)。