晶體銅微探針納米刻劃的分子動力學建模

張俊杰 孫 濤 閆永達 梁迎春 董 申

哈爾濱工業大學，哈爾濱，150001

0 引言

納米機械加工技術如金剛石刀具切削、金剛石磨粒加工以及金剛石微探針納米刻劃等的加工精度為納米量級，使加工形成的結構與器件的尺寸達到納米量級，是一種重要的由上而下的納米加工技術［1－2］。然而，缺乏對納米機械加工機理的理解極大地制約著納米機械加工技術的發展［3］。分子動力學（molecular dynamics，MD）模擬已被證實為研究納米機械加工機理的一個重要工具，國內外學者對納米機械加工的分子動力學模擬開展了大量的研究工作，并取得了很多有價值的研究結果［4－9］。然而，在他們的研究中，工件材料大多為單晶形式，對多晶材料納米機械加工的分子動力學模擬研究較少。多晶材料發生塑性變形時，晶界不僅會阻礙位錯運動，并且晶界自身也會發生變形，其變形機理遠比單晶材料復雜［10］。此外，前人研究主要依據宏觀統計量如切削力來分析分子動力學模擬的結果，而對工件材料內部缺陷演變的分析較少。缺陷的形核及運動是材料塑性變形的一種重要方式，并且工件材料的變形與宏觀統計量直接對應。因此，為了全面、深入地理解納米機械加工機理，需要使用分子動力學模擬并結合晶體缺陷分析技術來研究多種結構材料的納米機械加工機理。

納米機械加工的分子動力學模擬主要包含建立分子動力學模型、進行分子動力學模擬和處理分子動力學模擬數據。其中，建立包含具備基本物理特性的精確分子動力學模型決定著分子動力學模擬結果的正確性及精度。因此，本文針對上述不足之處，以晶體銅微探針納米刻劃為例，分別從分子動力學基本原理、晶體銅原子結構建模和晶體缺陷分析這三個方面開展了晶體銅微探針納米刻劃的分子動力學建模的研究。

1 分子動力學基本原理

圖1為晶體銅微探針納米刻劃的分子動力學模型（包含晶體銅工件和金剛石微探針）示意圖。晶體銅工件由具有面心立方（face centered cubic，FCC）晶格結構的銅原子組成，金剛石刀具由具有金剛石晶格結構的碳原子組成。工件在X方向、Z方向具有周期邊界條件，在Y方向具有自由邊界條件。同時，工件底部設置固定層原子以防止加工中工件的整體剛性移動。納米刻劃過程包含馳豫階段和加工階段。在加工階段中，微探針首先從工件表面上方以恒定的速度壓入工件內部直到指定的壓痕深度，然后沿著X軸負方向以恒定的速度進行刻劃直到指定的刻劃長度。

圖1 晶體銅微探針納米刻劃的分子動力學模型示意圖

在建立晶體銅納米機械加工的分子動力學模型時，需要從分子動力學基本原理角度出發設置經典牛頓運動方程的求解方法，描述模擬體系的勢能函數，模擬真實環境的系綜（ensemble）等參數。雖然分子動力學模擬已被證實為研究納米機械加工機理的一個有力工具，但是分子動力學模擬結果的正確性與精度受到分子動力學模型精度的影響。在本文建立的分子動力學模型中，采用velocity－Verlet法則來求解當前分子動力學模型中原子的位置、速度和加速度，采用 Mishin等［11］提出的EAM勢能函數的參數來描述晶體銅工件材料里的Cu－Cu原子相互作用，采用Morse勢能函數來描述晶體銅工件與金剛石微探針之間的Cu－C原子相互作用［12］。Morse勢能函數具有較大的結合能，導致工件表面與刀具之間存在一個很大的吸附力。在納米機械加工實驗中，由于工件表面很難做到原子級潔凈或者存在硫醇鈍化層，導致工件表面與刀具之間的吸附力顯著減小。因此，為了模擬鈍化層的存在對工件與刀具之間作用的影響，可以采用強排斥勢能函數來模擬剛性球形金剛石微探針，以消除吸附力對微探針與工件之間作用的影響。強排斥勢能函數的表達式為

式中，r為球形金剛石微探針中心與工件表面之間的距離；R為微探針的半徑；d為微探針的硬度［13］。

在晶體銅微探針納米刻劃的分子動力學模擬中，需要分別對馳豫階段和加工階段設定相應的系綜。在馳豫階段，采用等溫等壓（isothermalisobaric）系綜對模擬體系進行約束，要求模擬體系的溫度達到設定的溫度，系統壓力等于或接近于0。圖2所示為使用NPT系綜對單晶銅微探針納米刻劃的分子動力學模擬體系在目標溫度30K，目標壓力為0，下馳豫40×10－12s的過程中的溫度、壓力變化曲線。可以看到，在初期振蕩后，溫度和壓力均穩定在目標值，表明模擬體系馳豫完全。在隨后的加工階段中，由于工件材料的去除深度在數個或數十個原子層，因此只有很小一部分工件材料參與了加工過程，從宏觀上看，納米機械加工中工件材料體積的變化極小。在納米尺度下，熱量對材料變形的影響顯著增加，有必要研究納米機械加工中熱量的變化。因此，在加工階段可以使用微正則（microcanonical）系綜對模擬體系進行約束，要求模擬體系的體積和能量保持不變。

圖2 單晶銅微探針納米刻劃的分子動力學模擬體系的NPT馳豫過程

2 晶體銅建模

一個理想晶體是由晶胞在三維空間有序、周期排列而形成的。由于原子排列的周期性和疏密程度在不同的晶向上有顯著的差別，因此單晶銅的導電性、導熱性、彈性模量、屈服強度等物理性質均存在各向異性。納米機械加工是一個刀具與材料高度耦合的過程，單晶銅的各向異性直接影響工件材料的塑性變形及加工結果，因此可以通過創建具有不同晶面的單晶銅工件來研究單晶銅的各向異性對單晶銅微探針納米刻劃的影響。

多晶材料發生塑性變形時，晶界阻礙位錯運動，因此揭示位錯－晶界交互作用是研究多晶銅微探針納米刻劃加工機理的一個重要方面。為了排除晶界變形對工件塑性變形的影響，可以通過雙晶銅微探針納米刻劃的分子動力學模擬來揭示位錯－晶界交互作用。雙晶銅包含兩個具有不同晶體取向的晶粒以及連接它們的晶界。雙晶銅的原子模型可以根據重位點陣（coincidence site lattice，CSL）來創建：首先建立兩個相鄰的單晶銅晶粒，使它們具有相同的晶體取向。然后保持一個晶粒不動，讓另外一個晶粒繞著某一旋轉軸旋轉一定的度數。由于發生旋轉的晶粒的晶體取向不同于未發生旋轉的晶粒，因此兩個晶粒存在位相差，它們的晶格失配在晶粒之間形成晶界［14］。由CSL可知，不同的旋轉軸和/或旋轉度數得到的晶界結構也不同。按照相鄰晶粒之間的位相差的大小，晶界可以分為小角度晶界（相鄰晶粒的位相差小于10°）和大角度晶界（相鄰晶粒的位相差大于10°）。按照旋轉軸與晶界的關系，小角度晶界可以分為傾斜晶界（旋轉軸平行于晶界）和扭轉晶界（旋轉軸垂直于晶界）。還有一類特殊的孿晶界（twin boundary，TB），晶界兩側原子的對稱排列。圖3所示為根據CSL創建的包含不同晶界結構的雙晶銅工件的原子結構模型。

圖3 基于CSL創建的雙晶銅工件

多晶材料包含許多具有不同晶體取向的單晶晶粒，晶粒之間通過晶界連接。當晶粒尺寸進入到納米尺度后，由于納米晶體（nanocrystalline）材料的晶界體積分數（grain boundary volume fraction）顯著增大，因此納米晶體材料具有高強度、超塑性等獨特的機械性能［15］。納米尺度下，晶界對納米晶體材料在納米機械加工中的塑性變形及加工結果具有更為顯著的影響。因此，有必要創建納米晶體銅的原子結構模型，為開展納米晶體銅在納米機械加工的分子動力學模擬奠定基礎。目前，常用維諾圖（Voronoi diagram）來建立納米晶體材料的原子結構模型。設定P＝｛p1，p2，…，pn｝是一組隨機分布在一個平面或三維空間上的n個離散點，它將平面或三維空間分為n個Voronoi區域V（pi）。每個包含點pi的V（pi）是到pi距離最近點的集合，每個V（pi）對應一個晶粒［16］。然后在每個晶粒內根據指定或隨機的晶體取向填充FCC銅原子。不同晶體取向的晶粒存在晶格失配，在相鄰晶粒之間會形成晶界。當點pi分布在與笛卡兒坐標系某一坐標軸垂直的平面上時，得到的納米晶體銅具有柱狀（columnar）結構。圖4所示為基于維諾圖建立的納米晶體銅的原子結構模型，其中白色原子為晶界原子，灰色原子為FCC原子。

圖4 基于維諾圖建立的納米晶體銅的原子結構模型

3 晶體缺陷分析

受到晶體的生成條件、原子的熱運動、對晶體進行的加工過程、摻雜及其他條件的影響，實際晶體中往往存在偏離完整性點陣結構的晶體缺陷。晶體缺陷可以分為點缺陷（替位、間隙、空位）、線缺陷（位錯）和面缺陷（堆垛層錯和界面）等。晶體缺陷是隨著各種條件的改變而不斷演變的，如材料內部的位錯起源與增殖機制在外加應力作用下被激活，導致材料中的位錯密度隨著塑性變形的進行而增加。缺陷的形核及運動是金屬材料塑性變形的一種重要形式，通過晶體缺陷分析技術來精確辨別工件材料內部形成的晶體缺陷的種類和位置，是晶體銅納米機械加工的分子動力學模擬的一個關鍵之處。目前廣泛使用的晶體缺陷分析技術 有 徑 向 分 布 函 數 （radial distribution function，RDF）［17］、中心對稱參數（centro－symmetry parameter，CSP）［18］、共近鄰原子分析（common neighbor analysis，CNA）［19］、滑移向量分析（slip vector analysis，SVA）［20］、Ackland 鍵 角 分 布（bond angular distribution，BAD）分析［21］等。其中，共近鄰分析是一種廣泛應用于金屬材料的晶體缺陷分析技術。

在共近鄰分析中，假定有一個原子對包含原子α和β。對于一個給定的理想晶體，其結構可以用指數i、j、k、l組成的CNA圖表來辨別［19］：

（1）指數i。若原子α和原子β為近鄰原子，則i＝1；否則，i＝2。

（2）指數j。原子α和原子β共同享有的近鄰原子（共近鄰原子）的個數為j。

（3）指數k。原子α和原子β的共近鄰原子之間鍵的個數為k。

（4）指數l。當i、j、k相同時，原子α和原子β的共近鄰原子之間不同連接方式的個數為l。

如圖5所示，FCC結構里所有的共近鄰原子具有1421圖表形式。在密排六方堆積（hexagonal close packing，HCP）結構里，共近鄰原子具有1421和1422兩種圖表形式。體心立方（body centered cubic，BCC）結構的圖表分布不穩定，3/7的共近鄰原子具有1441圖表形式，4/7的共近鄰原子具有1661圖表形式。圖6為單晶銅擠壓后材料內部的缺陷結構的分子動力學模擬圖，原子根據計算得到的CNA值進行著色，其中，白色原子為處于表面或位錯等缺陷環境的Other原子，深灰色原子為處于堆垛層錯等缺陷環境HCP原子，灰色原子為處于完整晶格結構的FCC原子。為了更清楚地顯示缺陷結構，可以將未發生變形的FCC消除不顯示，如圖6b所示。可以看到，共近鄰分析可以清楚地辨別單晶銅工件中存在的缺陷結構，極大地促進了對晶體銅納米機械加工機理的深入理解。

圖5 晶體結構的CNA圖表示意圖［22］

圖6 單晶銅擠壓后工件材料內部的缺陷結構

4 結語

本文從分子動力學基本原理、晶體銅原子結構建模和晶體缺陷分析技術三方面介紹了晶體銅微探針納米刻劃的分子動力學建模。本文研究結果表明：通過合理地選擇經典牛頓運動方程的積分方法，精確描述體系的勢能函數，模擬真實環境的系綜，可以獲得到晶體銅納米機械加工的分子動力學模擬的精確原子模型。基于晶體結構、重位點陣和維諾圖等建立單晶銅、雙晶銅和納米晶體銅的原子結構模型，同時結合先進的晶體缺陷分析技術辨別納米機械加工過程中工件內部形成的缺陷的種類和位置，可以為系統地開展晶體銅納米機械加工的分子動力學模擬提供技術保障。

［1］Zhao X Z，Bhushan B.Material Removal Mechanisms of Single－crystal Silicon on Nanoscale and at Ultralow Loads［J］.Wear，1998，223（1/2）：66－78.

［2］Yan Y D，Sun T，Dong S，et al.Top－down Nanomechanical Machining of Three－ dimensional Nanostructures by Atomic Force Microscopy［J］.Small，2010，6（6）：724－728.

［3］陳日曜.金屬切削原理［M］.2版.北京：機械工業出版社，1993.

［4］Belak J，Lucca D A，Komanduri R，et al.Molecular Dynamics Simulation of the Chip Formation Process in Single Crystal Copper and Comparison with Experimental Data［C］//Proceedings of ASPE Annual Conference.New York，1991：100－109.

［5］Fang T H，Weng C I.Molecular Dynamics Simulations of Nano－lithography Process Using Atomic Force Microscopy［J］.Surface Science，2002（1/2）：138－147.

［6］Fang F Z，Wu H，Liu Y C.Modelling and Experimental Investigation on Nanometric Cutting of Monocrystalline Silicon［J］.International Journal of Machine Tools ＆ Manufacture，2005，45（15）：1681－1686.

［7］Mulliah D，Kenny S D，Smith R，et al.Molecular Dynamic Simulations of Nanoscratching of Silver（100）［J］.Nanotechnology，2004，15（3）：243－249.

［8］羅熙淳，梁迎春，董申.單晶鋁納米切削過程分子動力學模擬技術研究［J］.中國機械工程，2000，11（8）：860－863.

［9］郭曉光，郭東明，康仁科，等.單晶硅納米級磨削過程的理論研究［J］.中國機械工程，2008，19（23）：2847－2851.

［10］Jin M，Minor A M，Stach E A，et al.Direct Observation of Deformation－induced Grain Growth during the Nanoindentation of Ultrafine－grained Al at Room Temperature［J］.Acta Materialia，2004，52（18）：5381－5387.

［11］Mishin Y，Mehl M J，Papaconstantopoulos D A，et al.Structural Stability and Lattice Defects in Copper：Ab Initio，Tight－binding，and Embedded－Atom Calculations［J］.Physical Review B，2001，63（22）：224106.

［12］Zhu P Z，Hu Y Z，Ma T B，et al.Study of AFM－based Nanometric Cutting Process Using Molecular Dynamics［J］.Applied Surface Science，2010，256（13）：7160－7165.

［13］Ziegenhain G，Hartmaier A，Urbassek H M.Pair vs Many－body Potentials：Influence on Elastic and Plastic Behavior in Nanoindentation of FCC Metals［J］.Journal of the Mechanics and Physics of Solids，2009，57（9）：1514－1526.

［14］Tschopp M A，McDowell D L.Structures and Energies of 3Asymmetric Tilt Grain Boundaries in Copper and Aluminium［J］.Philosophical Magazine，2007，87（22）：3147－3173.

［15］Lu L，Sui M L，Lu K.Superplastic Extensibility of Nanocrystalline Copper at Room Temperature［J］.Science，2000，287（5457）：1463－1466.

［16］Du D S，Hwang F.Computing in Euclidean Geometry［M］.Singapore：World Scientific Publishing，1995.

［17］Matteoli E，Mansoori G A.A Simple Expression for Radial Distribution Functions of Pure Fluids and Mixtures［J］.Journal of Chemical Physics，1995，103（11）：4672－4677.

［18］Kelchner C L，Plimpton S J，Hamilton J C.Dislocation Nucleation and Defect Structure during Surface Indentation［J］.Physical Review B，1998，58：11085－11088.

［19］Honeycutt J D，Andersen H C.Molecular Dynamics Study of Melting and Freezing of Small Lennard－jones Clusters［J］.Journal of Physical Chemistry，1987，91（19）：4950－4963.

［20］Zimmerman J A，Kelchner C L，Klein P A，et al.Surface Step Effects on Nanoindentation［J］.Physical Review Letters，2001，87（16）：165507.

［21］Ackland G J，Jones A P.Applications of Local Crystal Structure Measures in Experiment and Simulation［J］.Physical Review B，2006，73：054104.

［22］Tsuzuki H，Branicio P S，Rino J P.Structural Characterization of Deformed Crystals by Analysis of Common Atomic Neighborhood［J］.Computer Physics Communications，2007，177（6）：518－523.