郭曉光 張 亮 金洙吉 郭東明

大連理工大學(xué)精密與特種加工教育部重點(diǎn)實(shí)驗(yàn)室，大連，116021

0 引言

目前在機(jī)械加工領(lǐng)域，分子動(dòng)力學(xué)仿真研究的對(duì)象幾乎都是無(wú)缺陷的理想晶體，沒(méi)有考慮材料實(shí)際存在的微小缺陷，也沒(méi)有涉及復(fù)雜的物理化學(xué)過(guò)程［1-5］。實(shí)際加工的材料即使經(jīng)過(guò)化學(xué)機(jī)械拋光，仍然或多或少存在一些空位、間隙原子、位錯(cuò)等缺陷，并且材料的性能在很大程度上取決于其內(nèi)部的缺陷，因此材料內(nèi)部缺陷對(duì)加工過(guò)程的影響不容忽視［6］。另外，材料內(nèi)部缺陷在超精密加工過(guò)程中受切削力、切削熱等的作用會(huì)逐漸發(fā)展變化，如微裂紋和微孔洞的形核、擴(kuò)展和匯合，從而導(dǎo)致材料的微觀(guān)結(jié)構(gòu)發(fā)生不可逆的熱力學(xué)耗散過(guò)程，造成材料表面、亞表面的損傷。因此，對(duì)考慮缺陷的單晶材料超精密加工過(guò)程進(jìn)行分子動(dòng)力學(xué)仿真將使分子動(dòng)力學(xué)方法更具有實(shí)際應(yīng)用意義。

本文建立了考慮空位缺陷的單晶硅磨削過(guò)程分子動(dòng)力學(xué)仿真模型，并進(jìn)行了仿真計(jì)算，分析了材料的變形、相變等，研究了超精密磨削加工機(jī)理中的一些重要問(wèn)題。

1 分子動(dòng)力學(xué)仿真

1.1 考慮空位缺陷的單晶硅磨削過(guò)程分子動(dòng)力學(xué)仿真模型

單晶硅磨削過(guò)程的三維分子動(dòng)力學(xué)模型如圖1所示。磨粒為金剛石原子，仿真中認(rèn)為磨粒是絕對(duì)剛性的，即不因外力作用而變形磨損。工件為單晶硅原子，仿真中，工件被分為牛頓層、恒溫層和固定邊界層［7］。牛頓層的原子運(yùn)動(dòng)可由牛頓方程來(lái)描述;恒溫層的引入是為了使磨削過(guò)程中產(chǎn)生的熱量及時(shí)傳導(dǎo)出去，通過(guò)不斷標(biāo)度原子速度保持該區(qū)域溫度恒定;固定邊界層的原子始終保持不動(dòng)，起到減小邊界效應(yīng)的作用，其內(nèi)部原子在仿真過(guò)程中不參與計(jì)算。

圖1 磨削過(guò)程分子動(dòng)力學(xué)仿真模型

隨機(jī)剔除原子來(lái)構(gòu)建相應(yīng)的帶點(diǎn)缺陷的單晶硅模型，在工件原子中剔除0.5%的原子，形成相應(yīng)的空位。存在的點(diǎn)缺陷會(huì)對(duì)晶體結(jié)構(gòu)產(chǎn)生一定的影響?；诘谝恍栽?，本文采用Material Studio中的CASTEP(Cambridge sequential total energy package)模塊來(lái)驗(yàn)證所構(gòu)建模型的可行性。CASTEP是Cerius2和Materials Studio的量子化學(xué)模塊之一，用密度泛函理論模擬固體、界面和表面的特性，是一種非常有效的研究點(diǎn)缺陷(空位、間隙、置換雜質(zhì))和擴(kuò)展缺陷(晶界和位錯(cuò))的性質(zhì)的模塊。圖2所示為無(wú)缺陷硅晶體的晶胞和存在一個(gè)點(diǎn)缺陷的硅晶體晶胞。理想晶體結(jié)構(gòu)硅晶體的參數(shù)為:晶胞棱長(zhǎng)a=b=c=0.543nm，夾角α=β=γ=90°。經(jīng)過(guò)CASTEP模塊的計(jì)算，存在一個(gè)點(diǎn)缺陷的晶胞的參數(shù)為:a=b=c=0.547nm，與理想晶格參數(shù)比較，偏差小于0.005nm;α=90.01254°，β = γ =89.98746°，偏差小于 0.005°?？芍?，當(dāng)點(diǎn)缺陷不大時(shí)，可以近似認(rèn)為晶體的結(jié)構(gòu)沒(méi)有發(fā)生變化，所以可以采用隨機(jī)剔除原子的方法構(gòu)建含點(diǎn)缺陷的單晶硅模型。

圖2 無(wú)缺陷和有空位缺陷的硅晶體晶胞結(jié)構(gòu)

1.2 勢(shì)函數(shù)的選取和力的計(jì)算

單晶硅是金剛石型晶體結(jié)構(gòu)，對(duì)硅原子勢(shì)能的計(jì)算應(yīng)考慮多原子價(jià)鍵之間的相互影響因素，故采用多原子系統(tǒng)的Tersoff勢(shì)函數(shù)u(rij)［8］對(duì)單晶Si－Si和金剛石C－C以及Si－C進(jìn)行計(jì)算，公式如下:

其中，u(rij)為原子i與原子j之間的勢(shì)函數(shù);fC(rij)為截?cái)嗪瘮?shù);fR(rij)為排斥項(xiàng)函數(shù);bij為低價(jià)函數(shù)，fA(rij)為吸引項(xiàng)函數(shù);rij為原子i與原子j之間的距離;θijk為原子i、j、k之間的夾角。其他參數(shù)的值見(jiàn)表1。

表1 Tersoff勢(shì)函數(shù)參數(shù)

勢(shì)函數(shù)確定以后，原子之間的作用力Fij就可以通過(guò)勢(shì)函數(shù)對(duì)rij求導(dǎo)得出，即

作用在第i個(gè)原子上的總原子力等于其周?chē)性訉?duì)該原子作用力的合力，即

1.3 系統(tǒng)的運(yùn)動(dòng)方程

由統(tǒng)計(jì)物理可知，對(duì)于一個(gè)由大量粒子組成的物理系統(tǒng)，其宏觀(guān)特性是這些粒子運(yùn)動(dòng)狀態(tài)的綜合反映，因此，分子動(dòng)力學(xué)仿真的核心問(wèn)題就是要計(jì)算所有粒子的運(yùn)動(dòng)規(guī)律和軌跡，故可假設(shè)系統(tǒng)中的每個(gè)原子都各自服從經(jīng)典牛頓力學(xué)規(guī)律，可直接用牛頓運(yùn)動(dòng)方程來(lái)描述:

式中，vi、ri分別為粒子i的運(yùn)動(dòng)距離和速度;t為時(shí)間;m為原子的質(zhì)量。

建立了上述模型后，根據(jù)現(xiàn)有計(jì)算機(jī)的能力，選擇適合的計(jì)算量，編制分子動(dòng)力學(xué)計(jì)算機(jī)仿真程序，仿真結(jié)果主要以三維形式顯示，工件尺寸為10nm×10nm×10nm，磨粒半徑為1.8nm，仿真溫度為293K，磨粒速度為20m/s，時(shí)間步長(zhǎng)為1fs(10－15s)，磨削深度為0.5430nm。

2 考慮空位缺陷的單晶硅磨削過(guò)程分子動(dòng)力學(xué)仿真結(jié)果分析

2.1 溫度分析

在模擬過(guò)程中，系統(tǒng)初始狀態(tài)設(shè)置為293K。模擬過(guò)程先對(duì)系統(tǒng)進(jìn)行弛豫，使系統(tǒng)達(dá)到平衡狀態(tài)，然后再進(jìn)行磨削過(guò)程的分子動(dòng)力學(xué)仿真［9］。圖3所示為磨削前進(jìn)行的16ps(1ps=10－12s)弛豫過(guò)程中的系統(tǒng)溫度變化曲線(xiàn)。初始狀態(tài)時(shí)，系統(tǒng)溫度急劇下降，在0.7ps左右時(shí)達(dá)到溫度最小值(約118K)。隨后系統(tǒng)溫度逐漸上升，至10ps左右達(dá)到相對(duì)穩(wěn)定值，隨后系統(tǒng)溫度在此附近振蕩，波動(dòng)幅度小于2K。隨著溫度的波動(dòng)減小，系統(tǒng)逐漸進(jìn)入平衡狀態(tài)。溫度初期的急劇變化一方面是由于初始模型是按照理想晶體結(jié)構(gòu)模型構(gòu)建的，存在一定的結(jié)構(gòu)不穩(wěn)定;另一方面由于工件中存在邊界層和自由表面，部分原子在弛豫過(guò)程中會(huì)重新排列，相互作用力較大。最終，系統(tǒng)通過(guò)恒溫層的溫度控制，在弛豫過(guò)程中逐漸達(dá)到平衡狀態(tài)。

圖3 弛豫過(guò)程中溫度的變化

磨削過(guò)程中系統(tǒng)牛頓層的溫度變化如圖4所示，磨削初始狀態(tài)時(shí)系統(tǒng)的溫度約為293K，隨著磨削過(guò)程的進(jìn)行，金剛石磨粒開(kāi)始逐漸接觸擠壓工件，原子間相互作用力增大，牛頓層的溫度急劇上升。隨著磨粒完全接觸擠壓工件，溫度上升趨勢(shì)變緩慢，逐漸趨于穩(wěn)定，最終達(dá)到穩(wěn)定切削溫度。這個(gè)過(guò)程可以理解為:當(dāng)磨粒開(kāi)始接觸工件后，磨粒與工件相互作用加劇，系統(tǒng)能量逐漸升高，原子間劇烈的相互作用使得牛頓層的溫度迅速升高;當(dāng)磨粒完全接觸工件后，恒溫層的調(diào)控作用使溫度上升趨勢(shì)變緩，通過(guò)調(diào)控作用，牛頓層溫度逐漸達(dá)到較為平穩(wěn)的切削溫度。在整個(gè)過(guò)程中，溫度出現(xiàn)波動(dòng)的原因來(lái)自于硅晶體中存在的少量點(diǎn)缺陷。在磨削過(guò)程中，新的原子進(jìn)入空位處，吸收大量能量，形成新的晶體結(jié)構(gòu)。恒溫層的溫度如圖5所示，在仿真中恒溫層的初始溫度為293K，恒溫層的溫度采取每10fs一標(biāo)定的方法，波動(dòng)范圍為283～303K，恒溫層用來(lái)調(diào)控系統(tǒng)的溫度。

圖4 牛頓層溫度

圖5 恒溫層溫度

2.2 勢(shì)能分析

圖6所示為弛豫過(guò)程中硅原子間勢(shì)能變化曲線(xiàn)，在弛豫過(guò)程中，系統(tǒng)的能量首先急劇增加，然后迅速降低，在恒溫層的調(diào)控下又逐漸增加，經(jīng)過(guò)不斷的調(diào)控，系統(tǒng)的能量趨于穩(wěn)定。這個(gè)現(xiàn)象可以解釋為:初始模型在結(jié)構(gòu)上存在大量空位點(diǎn)，這種不穩(wěn)定結(jié)構(gòu)會(huì)導(dǎo)致系統(tǒng)初期原子的劇烈振動(dòng)。隨著系統(tǒng)自身不斷的調(diào)控，部分原子重新排布，最終達(dá)到穩(wěn)定狀態(tài)。圖7所示為考慮空位缺陷的單晶硅納米級(jí)磨削加工過(guò)程的硅原子間勢(shì)能變化曲線(xiàn)，結(jié)果表明:隨著金剛石磨粒逐漸擠壓切削工件，磨削面積逐漸增加，產(chǎn)生塑性變形，導(dǎo)致硅原子間的勢(shì)能急劇增加，并伴有一定的往復(fù)波動(dòng);當(dāng)磨粒完全接觸工件后，系統(tǒng)的能量逐漸趨于穩(wěn)定，硅原子間勢(shì)能增加緩慢。磨粒與工件不斷相互作用，大量能量以晶格應(yīng)變能的形式貯存在單晶硅的晶格中，轉(zhuǎn)化為硅原子間的勢(shì)能，從而導(dǎo)致系統(tǒng)勢(shì)能逐漸增加。在磨削過(guò)程中，原子間勢(shì)能的反復(fù)波動(dòng)主要是由單晶硅晶格的變形、晶格重構(gòu)及非晶相變?cè)斐傻?。共價(jià)鍵晶體中位錯(cuò)的發(fā)生和運(yùn)動(dòng)要克服高能勢(shì)壘，會(huì)導(dǎo)致能量的增減，而圖7中勢(shì)能變化僅有很小的波動(dòng)，所以帶缺陷的單晶硅的納米級(jí)磨削加工過(guò)程中沒(méi)有出現(xiàn)明顯的位錯(cuò)。

圖6 弛豫過(guò)程勢(shì)能變化曲線(xiàn)

圖7 加工中勢(shì)能變化曲線(xiàn)

2.3 工件受力分析

磨削力是理解磨削過(guò)程的一個(gè)重要物理參數(shù)，它清晰反映了材料的去除過(guò)程。納米級(jí)磨削過(guò)程中磨削力的來(lái)源與宏觀(guān)磨削過(guò)程中磨削力的來(lái)源有很大區(qū)別:宏觀(guān)磨削過(guò)程中的磨削力是磨削過(guò)程中產(chǎn)生的切削力和摩擦力的總和;納米級(jí)磨削過(guò)程中磨削力的來(lái)源很簡(jiǎn)單，主要來(lái)源于單晶硅原子和金剛石磨粒原子之間的相互作用力。

圖8所示為考慮空位缺陷的單晶硅納米級(jí)磨削過(guò)程中，工件在X、Y、Z方向的受力曲線(xiàn)，從仿真結(jié)果看，X方向的力在波動(dòng)中逐漸增加，然后趨于穩(wěn)定，Y、Z方向的力僅在一定范圍內(nèi)波動(dòng)。通過(guò)力的變化分析可知，X方向的力在納米級(jí)磨削過(guò)程中起主要的切削作用，即磨粒的推擠作用在納米級(jí)磨削過(guò)程中占主導(dǎo)地位，材料的去除依賴(lài)于磨粒的推擠。在磨削過(guò)程中，工件原子晶陣在磨粒的推擠作用下發(fā)生原子鍵斷裂，原子晶格被破壞，形成非晶層。部分原子在磨粒Y方向力的作用下與已加工表層斷裂的原子鍵結(jié)合，重構(gòu)形成已加工表面的變質(zhì)層。處在磨粒前方的原子在磨粒的推擠作用下，被去除掉，形成切屑。從圖8可以看出，X、Y、Z方向的力都在往復(fù)波動(dòng)，有些地方波動(dòng)很劇烈。力的這種波動(dòng)與晶格的變形、晶格重構(gòu)和非晶相變有密切關(guān)系。當(dāng)磨削力增大并超過(guò)原子之間結(jié)合力臨界值且不足以形成位錯(cuò)時(shí)，原子點(diǎn)陣被破壞，導(dǎo)致原子鍵斷裂，成為非晶態(tài)原子。此時(shí)，磨削力會(huì)出現(xiàn)陡降至一個(gè)較低水平的波動(dòng)，而且這種作用力的波動(dòng)在整個(gè)磨削過(guò)程中不斷重復(fù)，因此可以認(rèn)為這是由非晶相變而產(chǎn)生的作用力波動(dòng)現(xiàn)象。

圖8 工件X、Y、Z方向受力曲線(xiàn)

2.4 磨削過(guò)程分析

磨削過(guò)程中瞬間原子位置如圖9所示。初始位置時(shí)，磨粒并未與工件接觸;5ps時(shí)，磨粒部分?jǐn)D壓工件，使工件受擠壓處原子重新排布;10ps時(shí)，磨粒已經(jīng)完全作用于工件，部分工件原子已經(jīng)離開(kāi)工件表面，堆積在磨粒前部。但由于磨削長(zhǎng)度較短，故并未堆積成切屑;15ps時(shí)，隨著加工面積增大，在刀具前方堆積成切屑，如圖9d所示。綜上所述，帶缺陷的單晶硅納米級(jí)的磨削過(guò)程可以解釋為:隨著磨粒逐漸向單晶硅片的擠壓，與金剛石磨粒接觸的最外層硅原子與金剛石原子間的作用力由引力轉(zhuǎn)化為斥力，同時(shí)外層原子又受到內(nèi)部硅原子的作用力。相比較排斥力，內(nèi)部作用力比斥力小，故斥力占主導(dǎo)地位。在磨粒的排斥作用下，磨粒下方的硅晶格發(fā)生了擠壓變形。隨著金剛石原子與硅原子之間距離的不斷減小，原子間的斥力增大，大量能量以晶格應(yīng)變能的形式貯存在單晶硅的晶格中，因而此能量也隨著力的增大而不斷增加。當(dāng)貯存在變形晶格中的應(yīng)變能超過(guò)一定值且不足以形成位錯(cuò)時(shí)，硅的原子鍵就會(huì)斷裂，規(guī)則的晶格結(jié)構(gòu)將被打破，原子排列逐漸變?yōu)闊o(wú)序狀態(tài)(非晶態(tài))，形成非晶層。另一部分原子離開(kāi)原位置，堆積在金剛石磨粒的前方形成切屑。

當(dāng)磨粒繼續(xù)向前加工逐漸離開(kāi)帶缺陷的單晶硅時(shí)，與金剛石磨粒接觸的最外層硅原子和金剛石原子間的作用力由斥力慢慢轉(zhuǎn)化為引力，同時(shí)它們又受到內(nèi)部硅原子的作用力，使非晶層的原子重構(gòu)，釋放了部分能量，從而達(dá)到新的平衡狀態(tài)。如圖9e所示，被擠壓下去的原子會(huì)出現(xiàn)一定的回彈，已加工表面并未形成一層平整的“磨痕”。在形成“磨痕”過(guò)程中，“磨痕”兩側(cè)也不會(huì)產(chǎn)生大量側(cè)向原子。這是由于硅是脆性材料，其脆性較大而斷裂韌性又較小，材料的彈性極限和強(qiáng)度相差不大。

圖9 磨削過(guò)程中瞬間原子位置圖

2.5 空位缺陷對(duì)加工的影響

硅晶體中存在的空位缺陷會(huì)對(duì)加工中系統(tǒng)溫度的變化和加工表面的質(zhì)量有一定的影響。圖10 所示分別為含0.5%、1.0%、1.5% 空位的缺陷硅晶體加工過(guò)程中系統(tǒng)溫度的變化曲線(xiàn)。加工初始的階段，三條曲線(xiàn)接近于重合，隨著加工的進(jìn)行，含空位越多的系統(tǒng)溫度變化越小。當(dāng)模擬到40ps時(shí)，系統(tǒng)溫度隨著空位數(shù)目的增加而逐漸降低。這是由于在模擬過(guò)程中，空位的存在會(huì)吸收大量的熱，空位越多，吸收的熱量越多，系統(tǒng)的溫度上升越慢。另一方面，空位的存在會(huì)在一定程度上影響到加工表面的質(zhì)量。含0.5%、1.0%、1.5%空位的晶體加工表面的回彈量分別為0.15nm、0.12nm和0.10nm。出現(xiàn)回彈是由于在磨粒推擠過(guò)程中，工件表層部分原子向下運(yùn)動(dòng)，形成非晶結(jié)構(gòu);失去磨粒作用后，此部分原子會(huì)受到原有晶格原子的作用，脫離晶格，在晶格上方重新形成已加工表面。在表層原子向下運(yùn)動(dòng)過(guò)程中，部分原子由于空位的存在，會(huì)進(jìn)入空位處，到達(dá)新的平衡狀態(tài)。當(dāng)失去磨粒作用后，原子的回彈量將會(huì)在一定程度減小。

圖10 不同空位數(shù)γ的系統(tǒng)溫度曲線(xiàn)

整個(gè)模擬加工過(guò)程中，空位缺陷的存在并未對(duì)加工過(guò)程中切削力的變化有較為明顯的影響。同時(shí)，在納米級(jí)仿真過(guò)程中，空位缺陷對(duì)加工表面質(zhì)量的影響不大。

3 結(jié)論

(1)基于第一性原理，建立并驗(yàn)證了含空位缺陷的單晶硅磨削過(guò)程的分子動(dòng)力學(xué)仿真模型。

(2)材料的去除過(guò)程可解釋為:在磨粒的推擠作用下，工件發(fā)生變形，去除原子堆積在磨粒的前方;當(dāng)貯存在變形晶格中的應(yīng)變能超過(guò)一定值時(shí)，材料表面的原子鍵斷裂，即完成了材料的去除。

(3)磨削過(guò)程中，磨粒前下方的硅晶格在磨粒的壓應(yīng)力作用下發(fā)生了變形，當(dāng)貯存在變形晶格中的應(yīng)變能超過(guò)一定值且不足以形成位錯(cuò)時(shí)，硅的原子鍵斷裂，晶格被打破，變?yōu)闊o(wú)序狀態(tài)，形成了非晶層。隨著磨粒不斷前移，損傷層逐漸向前、向深處擴(kuò)展，造成了單晶硅亞表面的損傷;同時(shí)，部分非晶層原子在壓應(yīng)力的作用下與已加工表層斷裂的原子鍵結(jié)合，重構(gòu)形成已加工表面變質(zhì)層。磨削過(guò)程中，無(wú)明顯的位錯(cuò)、裂紋等產(chǎn)生。

(4)磨削過(guò)程中，工件X、Y和Z方向的受力及硅原子間勢(shì)能都是往復(fù)波動(dòng)的，在空位集中處會(huì)出現(xiàn)較大波動(dòng)。力和勢(shì)能的這種波動(dòng)與晶格的變形、晶格重構(gòu)和非晶相變有密切關(guān)系。

(5)空位缺陷數(shù)目會(huì)影響到加工過(guò)程中溫度的變化，從而在一定程度上會(huì)影響連續(xù)加工的穩(wěn)定性。空位的存在會(huì)減小加工后表面的回彈，但對(duì)加工表面質(zhì)量的影響較小。

［1］Tanata H，Shimada S，Higuchi M，et al.Mechanism of Cutting Edge Chipping and Its Suppression in Diamond Turning of Copper［J］.Annals of the CIRP，2005，54(1):51-54.

［2］Komanduri R，Raff L M.A Review on the Molecular Dynamics Simulation of Machining at the Atomic Scale［J］.Proceeding of the Institution of Mechanical Engineering，Part B:Journal of Engineering Manufacture，2001，215(12):1639-1672.

［3］Tanata H，Shimada S，Anthony L.Requirements for Ductile－mold Machining Based on Deformation Analysis of Mono－crystalline Silicon by Molecular Dynamic Simulation［J］.Annals of the CIRP，2007，56(1):53-56.

［4］Han X S，Lin B，Yu S Y，et al.Investigation of Tool Geometry in Nanometric Cutting by Molecular Dynamics Dimulation［J］.Journal of Materials Processing Technology，2002，129(1/3):105-108.

［5］于思遠(yuǎn)，林濱，韓雪松，等.分子動(dòng)力學(xué)仿真技術(shù)在超精密加工領(lǐng)域中的應(yīng)用［J］.中國(guó)機(jī)械工程，2002，13(1):22-25.

Yu Siyuan，Lin Bin，Han Xuesong，et al.The Application of Molecular Dynamics Simulation to Ultraprecision Machining［J］.China Mechanical Engineering，2002，13(1):22-25.

［6］Shimada S，Tanaka H，Ikawa H.Atomistic Mechanism of Surface Generation in Micromachining of Monocrystalline Silicon［C］//Proceeding of the 1st International Euspen Conference.Bremen，1999:230-233.

［7］張偉文，郭鋼，黃云，等.納米磨削磨屑形成分子動(dòng)力學(xué)仿真研究［J］.中國(guó)機(jī)械工程，2011，22(2):127-131.

Zhang Weiwen，Guo Gang，Huang Yun，et al.A Molecular Dynamics Study of Chip Formation on Nanogrinding Processes［J］.China Mechanical Engineering，2011，22(2):127-131.

［8］Tersoff J.Modeling Solid － state Chemistry:Interatomic Potentials for Multicomponent System［J］.Physics Review B，1989，41:5566-5568.

［9］CheongW CD，Zhang L C.Molecular Dynamics Simulation of Phase Transformations in Silicon Monocrystals due to Nano － indentation［J］.Nanotechnology，2000，11(3):173-180.