以多尺度模擬與實驗驗證探討層狀結構的介面破壞行為(I)

(1)

行政院國家科學委員會專題研究計畫成果報告

以多尺度模擬與實驗驗證探討層狀結構的介面破壞行為(I)

計畫類別：個別型計畫計畫編號： NSC93-2211-E-002-048- 執行期間： 93 年 08 月 01 日至 94 年 07 月 31 日執行單位：國立臺灣大學土木工程學系暨研究所計畫主持人：陳俊杉計畫參與人員：黃仲偉、廖英博、賴家偉報告類型：精簡報告報告附件：出席國際會議研究心得報告及發表論文處理方式：本計畫可公開查詢

中華民國 94 年 10 月 30 日

(2)

中文摘要

本研究計畫為兩年期的研究，在第一年發展適用於奈米尺度的各項模擬工具與計算方法，在第二年與相關奈米實驗進行比對與探討。我們於本年度（第一年）成功的發展出兩套計算方法與工具，並以此探討差排對金屬材料力學行為之影響。計算方法與工具之一是適用於大尺度平行計算的分子動力模擬（molecular dynamics）、其二是透過自由度的適性選擇，耦合分子動力學與連體力學的跨尺度擬連體計算（quasicontinuum）。透過分子動力模擬與embedded-atom 古典勢能，我們探討在百萬顆（106）原子等級的元件尺寸於壓印實驗中所導致之差排生成與差排組合結構的微觀現象。在跨尺度擬連體計算方法中，我們成功的推導出內力計算所需的連結元素，並以一維與三維的案例，驗證此法的正確性。

一、前言

隨著電腦軟、硬體與計算速度的快速成長與實驗儀器解析度的大幅提升，微觀材料的計算與模擬已能逐漸與微觀的實驗觀察作詳細的比對與驗證。透過這些瞭解，我們可進一步的擁有預測與改善材料巨觀性質的能力，如材料的力學性質等（MRS 2000，MRS 2001）。材料的力學性質基本上取決於各種瑕疵（defects）在材料微結構中的交互作用，這些瑕疵包括差排、微裂縫、介面等。整體而言，本計畫的長期目標是希望能以計算模擬與實驗驗證，建構預測材料行為的模型與理論。本研究於本年度以原子尺度模擬探討差排對金屬材料力學行為的影響，此基礎的瞭解與釐清將有助於未來應用延伸至更複雜的材料與微結構系統，如多材料層狀結構、精密與智慧陶瓷、高分子材料等。此瞭解亦有助於未來應用延伸至多重物理現象，如力與光電磁熱耦合等。

二、研究目的

本研究於本年度發展適用於奈米尺度的各項模擬工具與計算方法，並以此探討差排對金屬材料力學行為的影響，研究的主軸以奈米壓印現象為主，奈米尺度模擬的工具包括發展創新的計算方法如平行分子動力計算與跨尺度擬連體法，以突破原子模擬的尺寸限制，達到預測與釐清差排對奈米元件的材料力學行為之影響。

三、以分子動力模擬探討奈米壓印行為

奈米元件的尺寸與厚度一般介於數奈米到數十奈米範圍，因此傳統的拉伸實

(3)

驗與硬度實驗已經無法適用，取而代之的是目前廣被使用的奈米壓印（nano-indentation）量測技術。奈米壓印量測技術起源於傳統的硬度實驗，但使用了更尖更小的壓印子（indenter）以產生奈米級的微小壓痕。藉由研究微小的壓痕及分析壓印時所產生的負載與壓印深度的關係，進而求出材料的彈性模數、硬度、韌性、強度等機械性質（Fischer-Cripps 2002）。過去幾年，以有限元素法為基礎的計算分析提供一有效方式探討壓印變形行為，Nix and Gao（1988）提出以應變梯度塑性理論（strain gradient plasticity）為基礎的材料組成律，成功的以幾何必要差排機制（geometrically necessary dislocations）預測與解釋介於數百奈米與數微米之間所觀察到的壓印尺寸效應。然而當尺寸小於100 nm 時，材料晶格瑕疵生成與互動行為顯得十分重要，使用以連體力學為基礎的有限元素分析，因為無法反應出晶格瑕疵對組成律的影響，故已無法正確模擬變形行為。欲進一步解決上述的問題，我們必須以原子為基本單位，模擬原子間的互動，由此取得系統內所有原子的位置及速度，並透過適當的後處理，進一步瞭解系統的性質。分子動力分析即是以此為主軸而發展，目前已被廣泛應用。分子動力分析能處理的原子數目受到電腦運算速度和儲存能力的限制。即使運用超級電腦進行平行處理，至今所能計算的原子數目只能達到 109 量級的原子數目（Abraham et al. 2002），這與現實物質動輒含有1023個原子的差異很大。針對此長度限制，我們進一步發展與跨尺度擬連體法，透過自由度的適性選擇，耦合分子動力學與連體力學。以下就分子動力分析的基本理論與成果作一簡述，跨尺度擬連體法的推導與驗證將於下一節說明。 3.1 基本理論 分子動力學一般用於小尺度（microscopic）系統的計算與模擬。該法首先需要建立一組原子的運動方程式，然後對系統中的每個原子透過數值積分求解取得系統內所有原子的位置及速度，即在相空間的運動軌跡（phase-space trajectories），再利用統計計算方法得到系統的靜態與動態特性，由此進一步瞭解系統的性質。考慮一含有N 個原子的系統，其對應之勢能函數為 V。則系統之拉

格蘭治運動方程式（Lagrangian equation of motion）可寫成下列形式： ( ) ( ) k k d L L dt q q ∂ ₋ ∂ ₌ ∂ ∂ 0 (3-1) 上式中qk為廣義座標，qk為廣義座標對時間的微分。L q q( , )k k 為拉格蘭治函數，假設第i 個原子在卡式座標下的位置與對應的動量分別為 ri與pi，則系統的拉格蘭治函數可以下式表達：

(4)

2 1 2 N i i i L m = =

∑

p −V i i (3-2) 上式中mi為第i 個原子的質量，而等號右邊第一項代表系統總動能，第二項則代 表系統總勢能。若將方程式（3-2）代回方程式（3-1），且令廣義座標 qk為各個原子的位置向量rk，則（3-1）式可改寫為： i i i mr =f (3-3a) i = ∇rL= −∇ V f r i i fi (3-3b) 上式中fi代表第i 個原子所受的作用力，其值等於勢能函數 V 對位置向量 ri作梯

度運算。對於卡式座標而言，系統第i 個原子的 Hamilton’s equation 可寫成（Allen

and Tildesley 1987）： / i = i m r p (3-4a) i = −∇rV = p (3-4b) 在三維卡式座標中，若要計算系統所含 N 個原子的運動軌跡，則可求解方 程式（3-3a）中 3N 個二階常微分方程組，或等價地求解方程式（3-4a、3-4b）中所含6N 個一階常微分方程組。在分子動力學的計算中為了兼顧計算效率與準確

度，一般採用有限差分法中的Verlet 或 Gere 演算法則來求解此微分方程組（Allen and Tildesley 1987）。

然而上述方程式只提到如何描述分子運動方程式，並未提到原子能量如何計算。在原子尺度下，原子的能量理應按量子力學計算。然而一般為簡潔起見，大都採用古典勢能（classical interatomic potential）。兩者的差異在於量子力學的勢能理論中原子間的能量是不可分割的，而古典勢能理論係利用近似或半經驗公式求得，原子間的能量被認為是可以分割（Curtin and Miller 2003）。對於不同的原子，由於鍵結形式並不相同，所以對應的古典勢能亦不相同。本計畫中所使用之古典勢能屬於Embedded-Atom Method (EAM)（Daw and Baskes 1984）簡述如下： EAM 能量計算除包括一般原子間的對勢能（pair potential）外，尚考慮圍繞在原子周圍電子所引致之嵌入能（embedded energy）。一般常用於金屬原子（如銅或鋁），第i 個原子所含的能量可寫成： 1 ( ) ( ) 2 i i ij ij ij j i E F ρ V r ≠ = +

∑

(3-5a) ( ) ij j ij j i r ρ ρ ≠ =

∑

(3-5b)

(5)

上式中F_i即為與電子密度有關之嵌入能，V 為第i 個原子與第 j 個原子間的對勢

能。而第i 個原子的電子密度

ij

ρ 其值為圍繞在其四周的原子電子密度的總和。本計畫探討的鋁金屬材料採用Ercolessi and Adams (1993) 所提出之鋁金屬之 EAM 勢能。在進行奈米壓印實驗中，除需利用EAM 勢能描述鋁奈米金屬材料中原子間的鍵結外，尚需要其他勢能來描述壓印子以及壓印子與鋁奈米金屬塊材間的作用力。一般為簡化計算起見，壓印子是假設為剛體（rigid body），原子之間的相對位置不隨壓印過程而改變，故可將壓印子簡化為一質點，僅需考慮此一質點與鋁奈米金屬塊材間的距離來判斷壓印子是否已接觸到鋁奈米金屬塊材。另一方面，壓印子與待測奈米元件間的作用力，文獻上至今仍無適當的描述方式。對於圓（球）形的壓印子，一般是假設壓印子質點與奈米元件中的原子之間的排斥能（repulsive potential）為： 2 ( _i) F = A R r− (3-6) 其中A 為常數，與壓印子的材質相關，R 與分別為壓印子半徑及球形壓印子中 心與元件中第i 顆原子間的距離（Lilleodden et al. 2003），若R 大於代表壓印子 與元件已經接觸。 i r i r 圖3.1 鋁奈米金屬塊材示意圖 3.2 壓印模擬成果與討論 如圖3.1 所示，本研究所模擬之鋁奈米塊材的尺寸和邊界條件如下：在 X 方向為 322Å，Y 方向 274Å，在 Z 方向 193Å。邊界條件在 X 方向和 Y 方向上為週期性邊界（periodical boundary condition），在 Z 方向底部為固定邊界（fixed boundary

Indenter Y 方向為 PBC Free Surface X 方向為 PBC X Y Z Fixed BC

(6)

condition）。壓印子所壓印的晶面分別為鋁金屬的(001)晶面和(111)晶面。使用的圓形壓印子的半徑大小為50Å，壓印子質點與奈米元件中的原子之間的排斥能由式(3-6)決定，使用的常數 A 值設定為 10。設置完畢後即可利用共軛坡降法進行能量最小化的疊代，求解平衡態之原子位置。在每一次疊代求解得到平衡態後，我們將壓印子向下移動0.1 Å 的距離再疊代求解平衡態，如此重複，一直到壓印子到達我們所指定的壓印深度值。 3.2.1 壓印鋁(001)晶面 壓印鋁奈米金屬塊材(001)晶面所得之壓印力與壓印深度曲線如圖 3.2 所示，整條曲線一開始呈現平滑的上升直到壓印深度到達 5.75Å(圖上 a 點處)時有一不連續的轉折點，文獻上認為這一不連續轉折點和差排缺陷(dislocation)的形成有關。差排缺陷形成時所伴隨的能量消逝反應到此曲線上，而a 點即是材料彈性行為的終點。為便於了解壓印模擬實驗進行時材料微結構的變化，可將模擬的原子依據特定之指標值作著色。在文獻上有多種指標值可供利用，而其中較廣為大家所使用的是Zimmerman et al. (2001) 所提出的滑動向量（slip vector），其定義 如下: 1 ( n s s x n ) α αβ αβ β α≠ = −

∑

− Χ (3-7)

其中，n 是最接近原子 α 的相鄰原子(nearest neighbor atoms)，ns為有滑動的原子

的數目，xαβ和Xαβ分別為在目前位置和參考位置上原子α 和原子 β 間的位置向量 差。參考位置指的是在材料不受任何應力作用時的各原子的位置。把各原子上滑動向量的長度當作指標值對原子作著色，可判斷那些原子有較大滑動值，即較有可能是缺陷區的原子，如此一來，材料的缺陷結構可以很容易的被辨別出來。表 3.1 分別將曲線上 a、b、c 和 d 所對應之材料微結構畫出，原子的顏色愈深(紅)，其滑動之程度愈大。在a 點之前並沒有任何差排缺陷產生，在 a 點時，曲線出現第一個不連續點，對應到材料差排缺陷的產生。隨著壓印深度的增加，差排缺陷逐漸變大變複雜。我們可觀察到差排缺陷全都產生在鋁金屬的{111}面上，這主要是因為{111}面是面心立方晶体（face center cubic lattice）的滑動面，與差排理論的預測相符合。此外我們也可觀察到在c 點的微結構中，有差排 cross slip 的

現象，即是差排從一滑動面滑移到另一滑動面的現象。觀察 c 點的微結構可發

現，有很明顯的二個面形成一個交角的情形，這就是差排corss slip 後留下的軌

跡。到達d 點後，差排缺陷的結構變得更複雜，可看的出有更多的差排由於壓印

(7)

圖 3.2 壓印鋁奈米塊材(001)晶面之壓印力與壓印深度曲線

表 3.1 壓印鋁(001)面之材料微結構

Displacement Top View Side View

(a) 5.75Å

(8)

(b) 7.15Å (c) 8.15Å (d) 8.85Å 3.2.2 壓印鋁(111)晶面 壓印鋁奈米金屬塊材(111)晶面所得之壓印力與壓印深度曲線如圖 3.3 所示，曲線一開始的表現和壓印(001)晶面的曲線一樣，平滑的上升直到壓印深度到達 4.25 Å (圖上 a 點處) 有一不連續點出現。觀察表 3.2 可發現這個不連續點也對應 到差排缺陷的產生，這現象和壓印(001)晶面時的現象類似。在 a 點，差排缺陷的產生使得曲線有一個明顯的load drop 產生，而之後曲線繼續上升，這是因為所產生的差排缺陷間的交互作用所產生的 dislocation locks 使得材料有硬化（hardening）的現象產生。曲線到 b 點又有一個不連續點，這對應到微結構的改變，差排缺陷後原本只在二個滑動面上活動，轉變成三個滑動面上都有差排活動

(9)

的痕跡。由表3.2 中可清楚看出，在 b 點時，這三個滑動面和材料表面交成一個四面體的結構。這三個滑動面間，兩兩交線各是一條差排線，在差排理論中這稱這種差排為Lomer dislocation，由於它是屬於不可移動的差排缺陷，所以會在兩兩滑動面間形成一個障礙，阻擋其它差排缺陷的移動，使的材料再次被硬化，曲線繼續上升。在c 點時，差排陷缺的結構和在 b 點時差不多，但由於所受的壓印力較b 點時大，所激起的滑動面明顯比 b 點時多。之後到達 d 點，不但所激起的滑動面比c 點時多，我們也觀察到，有其他的差排突破 Lomer dislocation 的阻擋，滑出四面體結構。圖 3.3 壓印鋁奈米塊材(111)晶面之壓印力與壓印深度曲線

(10)

表 3.2 壓印鋁(111)面之材料微結構

Displacement Top View Side View

(a) 4.25Å (b) 5.85Å (c) 7.55Å (d) 8.95Å

(11)

3.2.3 二晶面之壓印行為比較

在圖3.4 中，我們將壓印二個不同晶面之曲線畫出比較。小圖將二條曲線的

彈性部分畫出，我們可以看出壓印(111)晶面的曲線在壓印(001)晶面的曲線之上，這可以由晶面彈性異向性來說明。Kiely and Houston (1998)利用奈米壓印實驗分別比較和鋁同為面心立方晶體(FCC Lattice)的金，他們的實驗結果指出，金的(111)晶面的楊氏模數值大約比(001)晶面有高出 36%。我們推測同為面心立方晶體的鋁，晶面的彈性異向性之趨勢應該和金類似，(111)晶面的楊氏模數值高於 (001)晶面的楊氏模數值，這就解釋為何壓印(111)晶面的曲線在壓印(001)晶面的曲線之上。此外我們在圖上也可以看出壓印(111)晶面產生差排缺陷的時間點比壓印(001)晶面產生差排的時間點來的早。圖 3.4 壓印鋁奈米塊材(001)晶面和(111)晶面之壓印力與壓印深度曲線

四、改良式擬連體法探討

材料巨觀變形與破壞的現象肇因於微尺度的機制，因此若欲完整的瞭解與分析此效應，則不免必須耦合微、巨尺度之間的計算。針對奈米尺度的力學現象，一種可行的計算方式即是在小尺度中利用分子動力學描述原子的運動，用以觀察材料缺陷如差排等所造成的機制。而在大尺度中透過有限元素法來求解不規則幾何與不規則邊界的偏微分方程式，減少計算原子運動所需的自由度。

(12)

整合分子動力分析與有限元素分析是一兼具計算效率與準確度的多尺度計算方式，由於其應用相當廣泛，目前已有多位學者積極地從事這方面的研究 (Rudd and Broughton 1998; Broughton et al. 1999; Curtin and Miller 2003) 。然而此 方法在計算系統內力與能量上存在若干問題，無法正確計算有限元素區與原子區交界處之內力( Curtin and Miller 2003; Shenoy et al. 1999) 。在本計畫中，吾人提出過渡元素的觀念，利用過渡元素連接有限元素區與原子區，數值實例顯示透過此觀念可正確計算出系統的內力與能量。 4.1、基本理論 4.1.1 擬連體法理論 一般而言，分子動力學假設原子能量可由古典勢能（classical potential）計算，若系統包含N 個原子，且並未受到外力作用，則此系統的總勢能可表示為： ( ) Π u 1 ( ) _tot N _i i E = Π u = =

∑

E (4-1) 上式中u 是位移向量， 是此N 原子系統的總能量， 為第i 個原子的能量。 此系統的平衡條件即在於找到一組位移向量u 來最小化系統總勢能。 tot E E_i

基本上，擬連體法( Shenoy et al. 1999; Miller and Tadmor 2002; Knap and Ortiz 2001; Huang et al. 2004) 藉由 R 個代表原子（representative atoms）來近似 N 原子

系統，且R<<N。則系統總能量可改寫為： 1 R tot E n E_{α α} α= ≈

∑

(4-2) 上式中nα可視為『代表原子所代表的原子數』（或代表原子所對應之權重係數）。而重點即在於如何計算各個代表原子所對應之原子數nα以及準確的計算代表原子上的內力。首先，擬連體法將整個系統分為局部（local）與非局部（nonlocal）兩大區域。分類方式主要判斷其所處的變形場是否均勻（homogenous）而定。以圖 4.1 Lomer 差排為例，在差排核心部分由於變形較不均勻，視為非局部區域，該區域需要較多的代表原子。另一方面，在遠離差排核心之處被視為局部區域，由於變形均勻只需要少數的代表原子。該區域內其他被省略的原子則假設會依據有限元素的內差函數移動。如此可大幅減少計算所需的自由度。

(13)

圖4.1 Lomer Dislocation 示意圖 4.1.2 擬連體法理論能量計算 根據前述將整個系統分為局部（local）與非局部（nonlocal）兩大區域，則系統總勢能可進一步改寫為： 1 1 ( ) M NL tot E n Eα α α n Eβ β α= β= ≈

∑

F +

∑

(4-3) 上式中M 與 NL 分別為局部區域內所有元素的數目與非局部區域內所有代表原子 的數目，而 nα與 nβ分別為對應局部與非局部區域內代表原子能量的權重因子。由於非局部區域一般是以分子動力計算模擬，將每顆原子都視為代表原子，故nβ 可取為 1.0。另一方面，nα可直接計算每一元素內共包含多少原子，對於座落在元素邊界與頂點之原子，其權重則按比例計算 (Huang et al. 2004)。值得一提的是對每顆代表原子而言，若在其周遭Rcut距離（古典勢能作用範圍）之內可以觀察到其他代表原子，則這些代表原子都屬於非局部區域。若不滿足上述條件，則此代表原子屬於局部區域。另一方面，若某代表原子屬於非局部區域，但同時為某一元素之節點（意謂此元素為過渡元素，連接局部與非局部兩區域），則在計算能量時，將此代表原子歸於非局部區域內計算，而該元素計算權重時並不包含此一代表原子。對於局部區域內之代表原子（如遠離差排的區域），其能量Eα可依照Cauchy- Born 法則利用元素內變形梯度 Fα 近似計算( Shenoy et al. 1999)。而對於非局部區域之代表原子，其能量則考慮其周遭原子（Rcut之內的原子）的影響，按照分子動力學計算。根據前述的計算方式，由於不會有代表原子能量被重複計算，故可得到正確的能量。

若所考慮的系統屬於靜力且溫度為零的系統（static and zero temperature），

則系統的平衡條件即在於找到一組位移向量u 來最小化系統總勢能。故透過方程

(14)

最小值所對應之位移向量u，此即為系統平衡狀態下的位移。 4.1.3 擬連體法理論內力計算 在分子動力學中，若系統中任一顆原子移動，則在該原子周遭Rcut距離內的所有原子之間的交互勢能都會改變。既然系統總能量是個別原子能量的總和，則每顆原子所受到的作用力可簡單寫為： ,ij cut j tot i j r R i i E E < ∂ ∂ = − = − ∂

∑

∂ f x x (4-4) 上式中x i代表第i 個原子的位置向量，且實際計算時只需對該原子周遭的原子能 量進行微分與加總的工作。然而若依前節所述，將代表原子區分為局部（local）與非局部（nonlocal）兩大類計算能量的作法，在直接代入式（4-4）計算代表原子之間作用力時卻會產生問題：由於分子動力學與有限元素法計算能量之間的不諧和情形，在局部與非局部交界之間的過渡區會產生非物理現象的力（non-physical forces）( Curtin and Miller 2003; Shenoy et al. 1999; Miller and Tadmor 2002) 。本文則嘗試利用過 渡元素的觀念克服此一問題。首先仍將方程式（4-3）所得之近似總能量代入方程式（4-4）中，如此可得： 1 1 ( M NL i i i E E n α n β) α β α= β= ∂ ∂ ≈ − + ∂ ∂

∑

f x x (4-5) 上式等號右邊的兩項分別代表局部與非局部區域內代表原子上的作用力。對任一屬於非局部區域（且非過渡元素節點）的代表原子而言，其作用力計算方式與一般分子動力學並無差異（上式等號右邊第二項），可直接利用方程式（4-4）計算。而對於屬於局部區域之代表原子，其作用力可進一步化簡為： 0 1 1 0 0 1 ( ) M M i i i M E n n n α α α α α α α α α α α α α ε = = = ∂ ₌ ∂ ∂ ∂ _Ω ∂ ∂ ∂ ∂ = ∇

∑

F u x F u x P N X i Ω (4-6) 上式中εα為第α個元素的應變能密度，其對變形梯度Fα的微分為第α個元素的第一類 Piola-Kirchhoff 應力張量 Pα；而變形梯度 Fα對位移向量的微分 uα 則等於形狀函數之梯度∇ N （計算在 X₀ _α α之位置）；Ω0則為在Xα位置之單位原子體積。而應力張量Pα之計算方式可透過下式計算：

(15)

(

, 1 ij cut j j i j r R i E < ∂ ⎛ ⎞ = _⎜ ⊗ Ω

∑

_⎝ ∂ _⎠ σ x x x −

)

⎟ (4-7a) det( ) J_α = F (4-7b) _α 1 J α = α α− P F σ (4-7c) 上式中σ為 virial 應力張量，Ω為在 xα位置之單位原子體積，⊗為兩向量的張量積，Jα為變形梯度Fα的行列式值。然而對於同時連接局部與非局部區域之代表原子而言，由於兼具兩種特性（一方面是非局部區域的原子，另一方面是局部區域元素的節點），故不能單獨利用方程式（4-4）或（4-6）計算作用力。本研究將它視為過渡元素的節點，同時透過原子叢集（cluster）的概念( Knap and Ortiz 2001)，其作用力之計算方式如下： 1 , ( ) ik cut m i j k j k r R w = < ⎡ ⎤ ≈ _⎢ ⎣ ⎦

∑

f f N X_k _⎥ (4-7) 上式中 m 為與該代表原子相接之過渡元素總數，wj為每個相接元素所對應之權重因子。而中括號內的f k與N(X k)分別意謂著在此代表原子周遭 Rcut距離之內每顆原子的作用力（利用方程式（4-4）計算）與相對應之內插值（透過 X k計算形狀函數值）。此法可同時考慮到該原子的兩種特性。 4.2 擬連體法數值實例 4.2.1 一維問題 在本節中我們以一簡單的一維範例，說明本文所提出之方法不但能消除傳統擬連體法中於介面上所面臨的非物理力，亦能在高度非均勻變形場的條件下提供精準的分析結果。為了使非均勻變形的效果更明顯，我們設計了一個合成的對勢能，如圖 4.2 所示。比起傳統所用的對勢能其影響範圍更遠，Rcut的距離大約是三倍，能量函數採用如下型式： 12 2 1 6 4 ) ( r e e r = _r − _r + φ (4-8) 力的函數型式如下：

(16)

13 2 12 6 8 ) ( r e e r f = _r − _r + (4-9) 如圖4.3 所示，本題為一個類似傳統連體力學中集中載重的範例。首先布置一條長原子鏈，左端為原子區，右端為元素區，每十個原子選一個代表原子。將最左端的原子固定住，並於介面原子上施加一個向右側拉的集中力（point force）。 0 1 2 3 4 5 6 7 r( ) -2 -1 0 1 2 3 pot ent ia l/f o rc e potential energy force 圖4.2 一維遠距合成對勢能函數及力 50 61 50 61 圖4.3 一維範例示意圖在布置完畢後，即可開始進行疊代求解，求收斂後最小能量之平衡位置。在本例中我們也與傳統分子動力學的分析結果做比較，圖4 為位移的比較，從圖中我們可以看出本文提出之方法確實可以在較少的自由度下求出精度相仿的解。圖4.5 為交界面附近的應變分析結果，應變之計算方式如下： i i i i x u u ∆ − = −1 ε (4-10) 在圖4.5 中我們可發現，因長距勢能作用下所造成介面上應變的落差影響範圍較長。在與分子動力學分析結果的比較之下，足以證明本文提出之分析方法亦能在高度非均勻的變形區域內得到不錯的結果。

(17)

0 40 80 120 160 X( ) 0 0.2 0.4 0.6 0.8 displacement(Ux, ) MD QC 圖4.4 一維範例分析結果（位移） 50 60 70 80 X( ) -0.004 0 0.004 0.008 0.012 0.016 str a in MD QC 圖4.5 一維範例分析結果（應變） 4.2.2 三維問題 在本節中我們將說明三維 Lomer 差排分析之範例。本例中模擬的元素為鎳（Ni），其晶格間距（）為 3.52Å，模擬之差排形式示如圖 4.6，X 方向平行於差排滑動方向 0 a ] 1 10 [ ，Y 方向平行於滑動面法向量，Z 方向平行於差排線之方向 ] 010 [ ] 1 0 1 [ 。模擬尺寸於 X 及 Y 方向各延伸 200Å，Z 方向厚度為 4 倍之 2a₀，

(18)

如圖4.7 所示。 SLIP DIRECTION ] 1 10 [ DISLOCATIN LINE ] 1 0 1 [ SLIP PLANE ] 111 [ ] 100 [ ] 010 [ ] 001 [ X Y Z 圖4.6 三維 Lomer 差排範例示意圖圖4.7 也展示本例之代表原子選擇和生成網格之分布。差排核附近及滑動面兩區域因其屬於非均勻變形場，因此我們選擇較密之代表原子，亦即原子區範圍。越遠離這兩個區間則越接近均勻變形場，代表原子的選擇也越稀疏，亦即元素區範圍。本例中共選擇了8,529 個代表原子，其中 7,527 個為非局部代表原子， 1,002 個為局部代表原子。本例若以傳統之分子動力模擬需要 82,861 個原子，因此採用本文所提出之方法進行相同的模擬，可節省大約90%的自由度。接著利用高效率之Delaunay 網格生成軟體產生網格，在本例中共產生了 9,440 個元素。在設置完畢後即可利用共軛坡降法進行能量最小化的疊代，求解平衡態之分佈。如同上一節之範例，本例也與傳統分子動力學的分析結果相互比較，驗證其分析結果的優劣。圖4.8 及圖 4.9 為滑動面（[101]方向）上下之位移跳躍分析結果。圖4.8 為 X 方向（[101]方向）位移之跳躍，圖4.9 為 Y 方向（方向）位移之跳躍。 ] 010 [ 由本例之分析結果與分子動力學分析結果互相比較後，足以驗證本文提出之方法確實能於三維問題中利用較少之自由度有效地求解出精度相仿的結果。除此之外，於差排核附近屬於高度非均勻變形場，在本例中屬能量弛豫程度最大的區域，於Lomer 差排中將形成一個固定不動的差排核結構。如圖 4.9 所示，本文提出之方法確實能有效求解高度非均勻變形區域，與傳統分子動力學分析結果比較，差排核附近能量弛豫後所求得的解其精度亦在可接受之範圍內。

(19)

-10 0 10 [-1 0-1_]-100 -50 0 50 100 [10-1] -100 -50 0 50 100 [01 0] X Y Z 圖4.7 三維 Lomer 差排範例之代表原子分佈及生成網格示意圖 -80 -40 0 40 80 [10 1] 0 0.5 1 1.5 2 displace ment jum p( 2.5 ) MD QC 圖4.8 三維 Lomer 差排範例分析結果（[101]方向位移跳躍）

(20)

-20 -10 0 10 20 [10 1] -0.04 0 0.04 0.08 0.12 0.16 0.2 di splacem e nt jump( ) MD QC 圖4.9 三維 Lomer 差排範例分析結果（[010]方向位移跳躍）

五、研究成果自評

本研究計畫為兩年期的研究，在第一年發展適用於奈米尺度的各項模擬工具與計算方法，在第二年與相關奈米實驗進行比對與探討。我們於本年度（第一年）成功的發展出兩套計算方法與工具，並以此探討差排對金屬材料力學行為之影響。計算方法與工具之一是適用於大尺度平行計算的分子動力模擬、其二是透過自由度的適性選擇，耦合分子動力學與連體力學的跨尺度擬連體計算。透過分子動力模擬與 embedded-atom 古典勢能，我們探討在百萬顆（106）原子等級的元件尺寸於壓印實驗中所導致之差排生成與差排組合結構的微觀現象。在跨尺度擬連體計算方法中，我們成功的推導出內力計算所需的連結元素，並以一維與三維的案例，驗證此法的正確性。整體而言，本計畫的研究內容與原計畫相符，預期目標也超前原計畫第一年之目標。研究的成果也將於近期內整理，發表於國際學術期刊。另外也將積極與奈米實驗做進一步的比對探討，期能於不久的未來達到以計算模擬與實驗驗證，建構預測奈米材料行為的模型與理論。

參考文獻

Abraham, F. F., Walkup, R., Gao, H., Duchaineau, M., Rubia, T. D. D. L., and Seager, M. (2002). “Simulating materials failure by using up to one billion atoms and the world’s fastest computer: work-hardening,” PNAS, 99(9), pp. 5783-87.

Allen, M. P. and Tildsley, D. J. (1987). Computer Simulations of Liquids, Oxford

(21)

Broughton, J.Q., Abraham, F.F., Noam, B. and Efthimios, K. (1999), “Concurrent coupling of length scales: methodology and application,” Phys. Rev. B, 60, pp.

2391-2403.

Curtin, W.A. and Miller, R.E. (2003), “Atomistic/continuum coupling in computational materials science,” Modeling and Simulation in Materials Sciences and Engineering, 11, R33-R68.

Ercolessi, F. and Adams, J.B. (1994), “Interatomic Potential from First-Principles Calculations: the Force Matching Method,” Europhys. Lett. 26, pp.583-588

Fischer-Cripps A. C. (2002), Nanoindentation, Springer, New York.

Huang, C.-W., Liao, Y.-P. and Chen, C.-S. (2004), “An improved quasicontinuum method for coupling atomistic simulation and continuum at mesoscale,”

IUMRS-ICA 2004, Hsinchu, Taiwan, 6 pages, CD-ROM.

Kiely, J.D. and Houston, J.E. (1998), “Nanomechanical properties of Au (111), (001), (110) surface,” Phys. Rev. B, Vol. 57, NO. 19.

Knap, J. and Ortiz, M. (2001), “An analysis of the quasicontinuum method,” Journal of the Mechanics and Physics of Solids, 49, pp. 1899-1923.

Lilleodden, E. T., Zimmerman, J. A., Foiles, S. M., and Nix, W. D. (2003). “Atomistic simulations of elastic deformation and dislocation nucleation during nanoindentation,” Journal of the Mechanics and Physics of Solids, 51, pp.

901-20.

Miller, R. and Tadmor, E.B. (2002), “The quasicontinuum method: overview, applications and current directions,” Journal of Computer-Aided Materials Design,

9, pp. 203-239.

MRS Bulletin, Theme: Theory and Simulation of Fracture, May 2000.

MRS Bulletin, Theme: Computational Materials Science and Multiscale Modeling of Materials, March 2001.

Nix, W.D. and Gao, H. (1998), “Indentation size effects in crystalline materials: a law for strain gradient plasticity,” Journal of the Mechanics and Physics of Solids,

46(3), pp. 411-425.

Rudd, R.E. and Broughton, J.Q. (1998), “Coarse-grained molecular dynamics and the atomic limit of finite elements,” Phys. Rev. B, 58(12), R5893-5896.

(22)

Shenoy, V.B., Miller, R., Tadmor, E.B., Rodney, D., Phillips, R., and Ortiz, M. (1999), “An adaptive finite element approach to atomic scale mechanics -- the quasicontinuum method,” Journal of the Mechanics and Physics of Solids, 47(3),

611-641.

Zimmerman, J.A., Kelchner, C.L., Klein, P.A., Hamilton, J.C., Foiles, S.M. (2001), “Surface Step Effect on Nanoindentation,” Phys. Rev. Lett., 87, NO.165507.

以多尺度模擬與實驗驗證探討層狀結構的介面破壞行為(I)

行政院國家科學委員會專題研究計畫 成果報告