以分子動力學研究水分子在奈米尺度下的熱性質及噴印特性

國立交通大學

應用化學所

碩士論文

以分子動力學研究水分子

在奈米尺度下的熱性質及噴印特性

Studies of the molecular dynamics simulation

of water on the thermal property

and the performance of nanojet printing

姓 名： 郭 芝 穎

指導教授： 吳 建 興 教授

中文摘要

中文摘要

噴墨技術是一門重要的技術，且應用的領域相當廣泛，本研究應用分 子動力學模擬純水在奈米尺度下噴流現象。第一階段先以分子動力學模擬水 分子之物性，並且與文獻驗證。第二階段以金原子建立噴管系統，模擬水分 子在此系統中之噴出過程，以預測奈米尺度噴流之現象。 (一)第一階段以分子動力學模擬在不同密度下的水分子之狀態值，所採用的 密度分別為0.95g/cm3、0.975 g/cm3、1.0 g/cm3、1.025 g/cm3、1.05 g/cm3， 探討溫度與壓力的關係。 (二)在驗證水分子的物性之後，將水分子置入系統之中進行噴流模擬，探討 噴出過程中的流線分佈、壓力分佈、密度分佈等性質。並且探討在不同系統 溫度下對各項性質所產生的影響。在研究中發現系統中之水分子的蒸發會隨 著系統溫度增加而有變大的趨勢。 (三)最後，選擇最適合進行奈米噴流的參考溫度來做噴印實驗，經由對系統 的模擬分析，了解控制噴印變因。 關鍵字：分子動力學模擬、奈米尺度、奈米噴流。

Abstract

Abstract

T Thhee iinnjjeett iiss aann iimmppoorrttaanntt tteecchhnnoollooggyy aanndd hhaass bbeeeenn wwiiddeellyy aapppplliieedd iinn t thheeiinndduussttrryy..IInntthhiissrreesseeaarrcchh,,mmoolleeccuullaarrddyynnaammiiccssssiimmuullaattiioonn iissaaddoopptteedd ttoo s siimmuullaattee tthhee pprroocceessss ooff nnaannoojjeecctt.. OOuurr wwoorrkk ddiivviiddeedd iinnttoo ttwwoo ppaarrttss iiss aass f foolllloowwss.. ( (II)) MMoolleeccuullaarr ddyynnaammiiccss iiss aaddoopptteedd ttoo ssiimmuullaattee tthhee ggllaassss ttrraannssiittiioonn t teemmppeerraattuurree ooff ppuurree wwaatteerr aatt ddiiffffeerreenntt ddeennssiittyy,, aanndd tthhee ddeennssiittyy hhaass 0.95 g/cm3、0.975 g/cm3、1.0 g/cm3、1.025 g/cm3 aannd 1.05 g/cmd 3 ,, r reessppeeccttiivveellyy..TThhee ssiimmuullaatteeddppuurrppoosseeiissttooeexxpplloorreetthhee rreellaattiioonn bbeettwweeeenntthhee p prreessssuurreeaannddtthheetteemmppeerraattuurreeooffppuurreewwaatteerr.. ( (IIII)) WWee bbuuiilldd aa nnaannoojjeecctt ssyysstteemm ttoo pprreeddiicctt tthhee rreeoollooggiiccaall bbeehhaavviioorrooff p puurree wwaatteerr.. IInn nnaannoossccaallee,, tthhee bbeehhaavviioorr ooff ppuurree wwaatteerr iiss aa bblliinndd ssppoott.. BByy m meeaannss ooff aannaallyyzziinngg tthhee streamline ddiissttiirrbbuuttiioonn,, pprreessssuurree ddiissttaannccee d diissttrriibbuuttiioonn,, aanndd ddeennssiittyy ddiissttrriibbuuttiioonn,, wwee ccaann ddeeeeppllyy uunnddeerrssttaanndd tthhee b beehhaavviioorrooffppuurreewwaatteerriinnnnaannoossccaallee.. ( (IIIIII))FFiinnaallllyy,,ccoonnttrroollaannddssttuuddyynnaannoojjeettpprriinnttiinnggssyysstteemm.. K KeeyyWWoorrdd：：MMoolleeccuullaarrddyynnaammiiccssssiimmuullaattiioonn,,nnaannoossccaallee,,NNaannoojjeecctt..

目錄

目錄

中文摘要

... I

Abstract ... II

目錄

... III

圖目錄

... VI

表目錄

... XI

符號說明

... XII

第一章、緒論

... 1

1.1 研究目的與動機 ... 1

1.2 流體微觀理論 ... 8

1.3 分子動力學理論 ... 9

1.3.1 分子動力學方法介紹 ...9 1.3.2 分子動力學模擬之優點 ...10 1.3.3 分子動力學目前之應用 ...11

1-4 噴墨技術簡介 ...11

1-5 水勢能簡介... 16

第二章、噴流文獻回顧

... 20

第三章、研究方法

... 36

3.1 分子動力學理論 ... 36

3.1.1 分子動力學基本假設 ...36 3.1.2 分子動力學模擬流程架構 ...37 3.1.3 分子動力學系統初始化 ...38

目錄 3.1.4 分子動力學系統控制 ...40 3.1.5 分子動力學求解運動方程式 ...43 3.1.6 分子動力學簡化方法 ...45 3.1.7 分子動力學性質計算 ...52

3-2 勢能函數選擇 ... 53

3.2.1 分子間勢能: ...53 3.2.2 分子內勢能數據: ...55

第四章

模擬系統... 57

4.1 純水物性的模擬系統 ... 57

4.2 奈米噴流模擬系統... 59

4.3 奈米噴流系統中純水塊狀的製備 ... 61

第五章、研究結果與討論

... 64

5-1 水分子物性研究 ... 64

5.2 不同溫度下之噴流研究... 66

5.2.1 瞬時圖與現象觀察 ...67 5.2.2 液滴於板子上的圖形 ...69 5.2.3 速度分佈圖 ...70 5.2.4 平均壓力圖 ...72 5.2.5 密度分佈 ...72 5.2.6 流場圖 ...74 5.2.7 壓力瞬時分佈圖 ...75

5.3 噴印研究 ... 76

5.3.1 推出的水個數 ...77 5.3.2 瞬時圖 ...79 5.3.3 壓力分佈圖 ...81

目錄

5.3.4 分子個數分佈圖 ...82

第六章、結論與未來展望

... 84

參考文獻

... 85

圖目錄

圖目錄

圖

1.1 生物晶片之探針陣列及成品

[34]

... 3

圖

1.2 全彩 PLED 成品圖(2 種像素)

[35]

... 3

圖

1.3 錫球成長的專利噴嘴

[36]

... 4

圖

1.4 電子構裝中的錫球成長及成品

[37]

... 4

圖

1.5 液晶顯示器(LCD)

[38]

... 4

圖

1.6 TFT-LCD 的彩色濾光片

[39]

... 5

圖

1.7 噴流模擬

[42]

... 6

圖

1.8 縫合線模擬排向行為

[43]

... 6

圖

1.9 高分子流場分析模擬

[44]

... 7

圖

1.10 高分子壓印模擬

[45]

... 7

圖

1.11 高分子噴出模擬

[46]

... 8

圖

1.12 噴墨列印技術... 12

圖

1.13 連續式：位元偏向式... 13

圖

1.14 連續式：多重偏向式... 13

圖

1.15 熱泡式：頂噴型... 14

圖

1.16 熱泡式：側噴型

[2]

... 14

圖

1.17 壓電式：收縮管型... 15

圖目錄

圖

1.18 壓電式：彎曲型... 15

圖

1.19 壓電式：推擠型... 16

圖

1.20 壓電式：剪切型... 16

圖

1.21 水分子模型種類

[40]

... 17

圖

1.22 TIP4P 水分子模型示意圖... 18

圖

1.23 SPC/E 分子模型圖 ... 19

圖

2-1 不同頻率噴流情形 ... 21

圖

2-2 尾隨液滴往前端合併... 22

圖

2-3 液滴尾隨情況圖 ... 22

圖

2-4 液柱斷離過程圖（Koplik & Banavar） ... 24

圖

2-5 液柱斷離過程圖（Kawano） ... 25

圖

2-6 奈米噴流過程圖 ... 26

圖

2-7 液體噴流預測情形 ... 26

圖

2-8 雷射熔射法情況一 ... 27

圖

2-9 雷射熔射法情況二 ... 27

圖

2-10 奈米噴嘴陣列 ... 28

圖

2-11 奈米噴嘴的製造... 29

圖

2.12 剪切流場示意圖... 29

圖目錄

圖

2.13 收縮管流之模擬... 30

圖

2.14 相同孔徑不同溫度之噴流... 30

圖

2.15 相同溫度不同管徑之噴流... 31

圖

2.16 水分子與 DNA 的模擬 ... 32

圖

2.17 水分子結冰模擬

[41]

... 32

圖

2.18 奈米壓印製程模擬

[50]

... 33

圖

2.19 奈米擠出流場模擬

[47][51]

... 33

圖

2.20 奈米固態薄膜受剪切力碎裂之模擬

[49]

... 34

圖 2-21 模擬收縮擴張流埸

[52]

... 34

圖 2-22 模擬塊狀性質示意圖

[50]

... 35

圖

3.1 分子動力學模擬流程圖... 38

圖

3.2 Lennard-Jones 勢能方程式及作用力 ... 47

圖

3.3 鄰近列表法示意圖... 48

圖

3.4 Cell-linking 列表法示意圖... 49

圖

3.5 Morse 列表式勢能圖... 51

圖

4-1 系統示意圖 ... 59

圖

4.2 水分子初始位置(X-Y 平面的視角及側面圖) ... 61

圖

4.3 壓塊水分子(1.67 g/cm

3

) ... 61

圖目錄

圖

4.4 置於系統中平衡後成圓柱狀(密度約 1 g/cm

3

) ... 62

圖

4.5 塊狀放入系統... 62

圖

4.6 水分子自然膨脹於噴管中... 62

圖

5.1 水分子的 P-T 圖... 65

圖

5.2 文獻的 P-T 圖... 65

圖

5.3 比較 P-T 圖... 66

圖

5.4 瞬時圖... 68

圖

5.5 瞬時放大圖... 69

圖

5.6 液滴於板子上的圖形... 70

圖

5.7 噴嘴內部 X 方向速度分佈圖 ... 71

圖

5.8 噴嘴內部 Y 方向速度分佈圖 ... 71

圖

5.9 壓力趨勢圖... 72

圖

5.10 密度分佈圖... 74

圖

5.11 流場分佈圖... 74

圖

5.12 壓力瞬時分佈圖... 76

圖

5.13 管內與管外的原子數... 78

圖

5.14 噴頭形成液珠... 79

圖

5.15 瞬時圖... 80

圖目錄

圖

5.16 壓力分佈圖... 81

圖

5.17 水分子分佈圖... 82

圖

5.18 水分子個數與推進距離作圖 ... 83

表目錄

表目錄

表

1.1 TIP4P 水分子勢能參數... 18

表

1.2 SPC/E 水分子勢能參數 ... 19

表

3.1 三～五階 Gear Predict-Corrector Algorithm 參數表 .. 45

表

4.1 分子勢能參數表... 58

表

4.2 分子勢能函式表... 58

表

4.3 Au-H

2

O 間作用力參數... 60

表

4.3 噴流系統參數表... 63

表

5.1 不同溫度的系統平衡時間... 64

表

5.2 不同溫度的噴流系統... 66

表

5.3 不同量的噴印測試... 77

表

A.1 減縮單位換算表... 91

符號說明

符號說明

i r 粒子i 的位置 i v 粒子i 的速度 i a 粒子i 的加速度 k Boltzmann’s 常數 s k 鏈結彎曲常數 θ 鍵結彎曲角度 ε σ, L-J potential 參數 V 系統體積 N 系統粒子數 F 分子間作用力 P 壓力 R 壁面原子與晶格距離 U 勢能 rij 粒子i 與粒子 j 間的位置向量 rˆ 單位方向向量 rc 截斷距離 r0 鍵結平衡距離 T 溫度

第一章、緒論

1.1 研究目的與動機

放眼望去，奈米已成為提高產品身價的重要形容詞，似乎只要是 奈米級的各種產品，就會有好品質，當然價錢也會變高，重要的是， 消費者依然願意花較高的費用換取較好的品質，例如: 奈米光粉底就 比一般粉底貴一至二倍，卻廣受女性消費者的喜愛，其受歡迎的原因 是奈米光粉底的粒子十分微小，可以將粗大毛孔完全填滿，讓臉蛋呈 現光滑透徹，這是就經濟而言 ; 對於醫療上的研究，也有極大的貢獻 跟期待，例如: 利用奈米載體的投藥標識，期望在未來的日子裡，藥 物可以只均勻分布於患部，而不是靠服藥的方式，造成高劑量藥物均 勻分佈於全身，可將服藥的副作用降到最低 ; 當然，也有運用在安全 上的實例，如: 奈米防爆隔熱紙，因為能與玻璃完全緊密地接合，可 使駕駛人發生車禍時，不被碎玻璃噴刺到，將傷害降到最低 ; 奈米科 技的開發與應用，確實替現在的科技產業帶來無限商機。 還有工業上的研究，隨著精密加工技術的不斷進步，微小化、 輕量化已成為趨勢。在此趨勢下奈米科技倍受矚目，因為奈米科技符 合現今工業上的自動化、微小化、降低成本、降低時程等要求。因此 全世界幾乎都在奈米科技上展開相關研究，也因為其創造性與影響力 十分廣泛，奈米科技成為 21 世紀的革命性產業，將使人類生活發生 劇變的科技技術，不僅在現有的科學與技術領域創造出新事物的可能 性，更打破過往的應用限制。 當全世界將大量資金投入奈米科技研究的同時，奈米尺度物質的 製備與量測就成為最重要的關鍵，但人類對於奈米科技還在摸索中，

無論是製作或測量上都得花費大量的時間和金錢，對於奈米科技研發 過程中會浪費許多資源。所以期望運用電腦模擬使奈米科技研發過程 變的有效率且節省資金。 當電腦模擬技術的概念逐漸形成，就會被廣泛使用在各個領域。 近年來常見的模擬有實體模擬、流場模擬、噴嘴模擬、模流分析模擬 等，這些模擬技術都是為了降低研發成本。 其中,噴墨列印技術能夠以低成本、快速、高精確的定量微液滴體 積，並且透過電腦協助，將液滴精確的置於定位，是一個極具潛力的 材料微分配技術，本研究將藉由分子模擬的方式，進行乙二醇水溶液 的奈米噴流時液滴形成及其穩定性現象探討，藉此來有效控制液滴大 小及位置。 奈米噴流行為特殊的現象使得預測變得困難，將利用分子動力學 (Molecular Dynamics, MD)理論，以分子為基本元素做模擬預測，以求 得更精準的奈米尺度噴流行為。 隨著這些科技的發展微小化將是其發展趨勢，為了面對更細微的 奈米尺度加工技術，本研究將採用分子動力學之方法來分析奈米尺度 下的噴流行為。微流體的應用相當廣泛，目前較受矚目的技術是微射 出成型技術(Micro molding)、生物晶片(Bio-chip)等精密高科技產業， 當然也包含本研究的微噴流技術(Nanojet)。 微噴流技術:另一微奈米尺度最具發展潛力的技術是微噴流技 術，因其具有低成本、快速、高精確微定量等絕佳的量產化條件，且 隨著所能控制的精度增加，其應用層面也相對不斷增加，此技術除了 在液滴控制、定位控制上的應用外，並且提供了一個極具潛力的材料

微分配技術，從早期的噴墨印表機、IC 封裝的應用到現在的生物晶片 (Bio-Chip)、有機薄膜電晶體(Organic Thin-Film Transistor, OTFT)、高 分子電激發光二極體(Polymer Light Emitting Diode, PLED)等方面 上，都一再的顯示此技術的優越性。 噴墨技術潛在應用面相當的廣泛，在此簡介以下應用:生物晶片之 引子和探針陣列、全彩PLED 製程、電子構裝中的錫球成長、TFT-LCD 的彩色濾光片。 圖1.1 生物晶片之探針陣列及成品[34] 圖 1.2 全彩 PLED 成品圖(2 種像素) [35]

圖1.3 錫球成長的專利噴嘴[36]

圖1.4 電子構裝中的錫球成長及成品[37]

圖1.6 TFT-LCD 的彩色濾光片[39] 然而，分子動力學最大的限制在於電腦的運算能力，當處理的系 統分子數目過於龐大時，所需的模擬時間將會過久，而浪費大量的時 間等待結果，不符何其模擬的初衷，因此過往一直未能擴大分子動力 學的應用領域。而近十年來因個人電腦效能快速提升、記憶體容量一 再擴充，使得分子動力模擬奈米尺度行為已成為可行之範圍。 國立清華大學的 CAE 實驗室，張榮語教授近年來所指導的學生 全力研究分子動力學，成果如圖1.7~圖 1.22，帶領國內利用分子動力 學進行模擬的研究。

圖1.7 噴流模擬[42]

圖 1.9 高分子流場分析模擬[44]

(1) (2) (3) (4) (5) (6) (7) (8) (1) (2) (3) (4) (5) (6) (7) (8) 圖 1.11 高分子噴出模擬[46]

1.2 流體微觀理論

本研究模擬奈米尺度噴流行為中的流體，選擇的流體是純水，所 以先分析了解流體於模擬時，分為哪些類型，主要被分為兩大類型， 分別為： 1. 分子流場模型(微觀):以分子為基本單位，根據分子間的交互作 用及能量關係來分析流體之行為，藉由直接測量和統計平均 來求取流場之性質。 2. 連續流場模型(巨觀):以質量、能量、動量等守恆定律發展出的

連續方程式，如Navier-Stokes 方程式、Euler 方程式等，將流 場視為連續體作計算。

因為巨觀時，做一些忽略的假設，而導致模擬現象不真實，而分子 流場模型就能較真實的描述流場行為，總言而之，微觀流體行為有別

於 往 常 巨 觀 下 的 連 續 體 ， 以 往 忽 略的 因 素 ： 如 表 面 張 力(Surface

Tension)、凡得瓦爾力(Van Der Waals Forces)、靜電力(Electrostatic Force)、空間力(Steric Force)，與介面之間的現象如邊界滑動…等，將 重新被考慮其效應。 所以初步研究，必須先將純水溶液在奈米尺度下的物性及流動行 為，利用分子動力學(MD)理論分析透徹，所以接下來先介紹分子動力 學理論。

1.3 分子動力學理論

當研究與應用的尺度隨著科技的進步而逐漸縮小時，過去適用於 巨觀模擬分析的連續方程式，因為其忽略許多微小尺度所存在的效 應，使得適用範圍將逐漸不足以涵蓋現今發展的奈米科技技術，因 此，以分子觀點出發的分子行為模擬就成為一項重要的分析方法。

1.3.1 分子動力學方法介紹

分子模擬的技術可以區分為兩種，一則是以統計隨機處理預測的 蒙地卡羅(Monte Carlo, MC)，一則是根據牛頓運動定律預測的分子動 力(Molecular Dynamics, MD)，也因為本質方法上的差異使其所適用 的物理行為也不盡相同。

蒙地卡羅分子模擬(MC)方法是採用隨機亂數的方式，因此無法真 實模擬物質的動態現象，大多數被使用在計算物質靜態平衡或穩定狀 態的性質探討，相較於分子動力模擬具有較快速找到穩態結構的優 點；相反的，分子動力模擬所需耗費的時間雖然較長，但是可以真實 模擬分子的動態行為，提供流動分子的完整行為描述，藉此可用來了 解流體在奈米尺度下的流動行為。 分子動力學的理論是將分子視為一個質點，藉由計算分子間受力 情形進而模擬出分子運動軌跡的科學，可用來模擬真實狀態的分子行 為並且藉此來求得系統的各種特殊性質；分子間的作用力計算方式採 用古典力學(牛頓運動定律)的物理行為。 分子動子學最早被提出來約在 50 年代，由於分子動力學模擬需 要龐大的計算量，而當時提出此理論的年代背景中，並不具有足夠運 算速率的電腦來進行此研究，僅能模擬簡單而少量的系統行為，近十 年來由於電腦科技快速蓬勃的發展，其處理資料量和運算速率遠遠勝 於從前，分子動力學技術因而又逐漸受到重視成為許多科學家的焦 點。

1.3.2 分子動力學模擬之優點

分子動力學最大的優點在於模擬的單元為奈米級尺寸的分子，可 精準的探討分子間交互作用的行為。然而研究實驗所能控制的精度還 未達分子級的程度，對於此極微小的研究仍有很大的困難。 即便觀測實驗成果，雖然目前已有不少測量的儀器設備，如:原子

力 顯 微 鏡(Atomic Force Microscope ,AFM) 、 掃 描 式 電 子 顯 微 鏡

(Scanning Electron Microscopy ,SEM) 或 穿 透 式 電 子 顯 微 鏡 (Transmission Electron Microscopy ,TEM)等，但除了儀器價格昂貴(百

萬~千萬)外，在量測上也有相當多的限制，且仍無法觀測奈米尺度的 動態行為。 此外奈米材料的製作和控制仍有許多困難，假使能利用分子動力 學進行模擬分析，不單系統參數控制簡單，就連現象觀測也都能藉由 電腦繪圖而容易了解，除此之外，對於某些超乎常態的高溫、高壓也 都能輕易的進行假設分析探討，除了省時、省力外，也能有效節省實 驗的成本。

1.3.3 分子動力學目前之應用

分子動力學模擬在除了被使用在探討分子、原子的材料性質或結 構物性外，在近年來也開始被使用在微小特殊元件的探討研究上。 在目前分子動力學使用的最多且最成熟的，是在於分子生物上的 探討，例如:蛋白質結構分析、DNA 基因片段、生物酵素反應機制等， 可以藉由分子動力學模擬瞭解其結構形成與反應物性，藉此來研究反 應機制與探討反應條件之影響。再加上近年來熱門的奈米材料也可藉 由分子動力學模擬來預測其形成之結構與條件，藉以幫助一些現象的 探討，如:奈米碳管、碳球(C60)之形成、微奈米流道行為或微結構電 子顯微鏡等，不只使得微奈米尺度的研究更便於瞭解，更可提供一個 具體的行為過程來進行分析。在過去十年來，分子動力學也經常被使 用在許多巨觀現象的研究，諸如:破壞力學、相變化行為、擴散現象、 滲透現象、雷射激發，或一些材料性質、熱力參數的探討，可以幫助 現象的細部結構研究，並可透過模擬來找尋適當的參數條件。

1-4 噴墨技術簡介

噴墨列印技術依照所使用的技術不同，其設計也相當多樣化，且

都有相當廣泛的應用潛力，圖 1-12 即為依不同技術所做的分類。一 般來說主要可分為連續式（Continuous）和供需式（Drop-On Demand, DOD）兩種方式。 圖 1.12 噴墨列印技術 連續式主要是利用墨水源連續產生墨滴，接著經過電場使得墨滴 帶有電荷，再藉由靜電或電磁方式使得墨滴上下偏向，將所要使用的

墨滴偏向列印的基材上，而不用的墨滴則偏向集墨管回收；依照偏向 方法的不同，可以分為位元偏向式（Binary Deflection）和多重偏向 式（Multiple Deflection），分別如圖1-13 和圖 1-14 所示。其中兩者的 優點皆為列印速度快且適用的化學藥品多；缺點為初始系統貴、解析 度低，且只適用低黏度墨水，這一類的應用主要是戶外的大型海報。 圖1.13 連續式：位元偏向式 圖1.14 連續式：多重偏向式 目 前 最 廣 泛 使 用 的 噴 墨 列 印 技 術 為 供 需 式 ， 而 其 中 熱 泡 式

（Thermal Bubble）和壓電式（Piezoelectric）為市場上的主流。熱泡 式是利用加熱噴嘴內部的薄膜電阻而產生高溫，使得墨水瞬間沸騰形 成氣泡，再藉由壓力上升產生的力量，將墨水由噴嘴噴出；熱泡式之 優點為設計簡單、價格便宜及速度快，缺點則是受熱應力的作用，噴 墨頭容易損耗。熱泡式主要可分為頂噴型（Roof-Shooter）和側噴型 （Side-Shooter），如圖 1-15 及圖 1-16 所示。 圖1.15 熱泡式：頂噴型 圖1.16 熱泡式：側噴型[2]

壓電式噴墨技術是利用壓電材料（Piezoelectric Material）施加電 壓後會產生形變，進而推擠噴嘴內墨水產生高壓而將墨滴噴出；優點 為無熱應力反覆作用因此不易損壞、墨滴大小容易控制，缺點為製造 技術困難且昂貴。依不同壓電產生的形變機制，壓電式又可分為收縮

管型（Squeeze Tube Mode）、彎曲型（Bend Mode）、推擠型（Push Mode）

及剪切型（Shear Mode），分別如圖 1-17、圖 1-18、圖 1-19 及圖 1-20 所示。

圖 1.17 壓電式：收縮管型

圖1.19 壓電式：推擠型 圖1.20 壓電式：剪切型

1-5 水勢能簡介

介紹水分子勢能模型，水分子為最常使用的溶劑，因此也發展的最 早，目前使用在水分子的勢能函式相當繁多，依水分子模型種類如圖 1.21，可以被區分為四大類，分別為：

a.五點水分子勢能(5-site Water Model)

b.四點水分子勢能(4-site Water Model) TIP4P potential(1983)

c.三點水分子勢能(3-site Water Model)

SPC potential(1981)、SPC/E potential (1987) d.Fitting to ab initio Water Model

MCY potential (1976)、CC potential (1984)

圖1.21 水分子模型種類[40] 由於水分子勢能模型種類相當繁多，發展上也各有優缺點，在過 去文獻中有許多水勢能的相關比較，一般而言，SPC、TIP4P 水勢能 模型是較經常被使用的水勢能模型，由於這兩種模擬方法較為容易， 而模擬出的相關性質也都與實際數值相接近。 一、TIP4P potential

TIP4P potential為The transferable Intermolecular Potential function

圖1.22 TIP4P 水分子模型示意圖 除了氫原子和氧原子外，在質心部位增加了一個帶電點帶有負電 荷，並且將氧原子的電荷省略，最重要的假設為此水分子勢能模型將 氫氧的鍵距以及鍵角都加以固定，故此水分子並不會有鍵角的變化和 鍵長的拉伸，以如此來代表整個水分子的勢能行為，也因此有計算上 較快速的優點，其勢能函式如下所示： ) ( ) ( 2 6 12 − +

∑∑

= µ µ µ v v v ij ij ij ij ij r e q q r C r A r U (1– 1) 所使用的參數如表1.1。 表1.1 TIP4P 水分子勢能參數 TIP4P 參數值 rOH(nm) 0.09572 rOM(nm) 0.015 ∠ HOH(o_{) 104.52} Aij(kcal/mol)Å12 6 x 105 Cij(kcal/mol)Å 6 610 qH(C) 0.52e qM(C) -1.04e

二、SPC/E potential 為三點式水分子勢能之一，由Berendsen等人於1981年所提出來的 勢能函式，所使用的水分子勢能模型是以三個原子的點位置來進行模 擬，其勢能形狀如圖1.23。 圖 1.23 SPC/E 分子模型圖 三個點分別為帶正電的兩個氫原子和帶負電的氧原子，其勢能函 式如(1-2)所示： ) ( ) ( ) ( 4 ) ( 2 6 12 ₊

∑∑

⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ − = µ µ µ σ σ ε v v v ij ij ij r e q q r r r U (1 – 2) 與先前的勢能函式是相同的，只在於分子的帶電量以及每個分子 相互間的位置和參數的改變，其參數如表1.2。 表 1.2 SPC/E 水分子勢能參數 SPC/E 參數值 rOH(nm) 0.100 ∠ HOH(o_{) 109.47} σ (nm) 0.3166 21 10− × ε (J) 1.0797 qH(C) 0.4238 e qO(C) -0.8476 e

第二章、噴流文獻回顧

液體噴流（Liquid Jet）的研究已有相當的歷史，並且在很早就已 經證實可藉由適當的控制，讓噴流斷裂（Breakup）產生液滴。早期

由Lord Rayleigh 提出線性理論（Linear Theory），他假設一無限長，

起始為靜止，不考慮黏度和周圍氣體影響，且為不可壓縮的液體噴 流，當給予一軸對稱的擾動後，半徑隨著表面張力(γ)之擾動的變化 為

(

t ikz

)

r =1+

κ

exp

α

− (2-1) 而擾動的增長速率(α)為

( )(

)

( )

2 1 0 2 1 2 1 3 1 ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ − ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ = k I k k k I a

ρ

γ

α

_(2-2) 其中κ為所給予的擾動，α為擾動的增長速率，k 為波數，γ為表 面張力，ρ為流體密度，a 為液柱的初始半徑。 其中長度和時間由初始半徑和 a/vc 將其無因次化，vc 為毛細波 （Capillary-Wave）速度，

(

)

12 a v_C = γ ρ ，由上可知，整個擾動是穩定或 非穩定主要根據波數k＜1 或 k＞1 而定，或者說λ 2a >π 或λ 2a<π （λ為波長）而定，而當波數等於0.697 時，此時α為最大值；若波 數小於 1，則由上述的分析可知，當

κ

exp

( )

α

t =1時噴流將會斷裂。 而 Rayleigh[1]在之後的研究發現，有時在主液滴和噴流間會有尾隨液 滴（Satellite Drop）出現，並用線性理論去解釋，不過後來証明這種 現象是屬於非線性的行為。

液滴生成現象相當的複雜，且液滴形成有時會伴隨著尾隨液滴的 產生，因此開始有學者針對這種現象進行研究，所以有非線性的探討 出現；例如 1968 年 Yuen[2]和1970 年 Nayfeh[3]分別藉由理論推導來研 究非線的影響；而在1970 年 Goedde 和 Yuen[4]由實驗來觀察非線性的 影響，他們發現尾隨液滴會隨著波數的減少而增大,由圖中可看出。 圖2-1 不同頻率噴流情形 1971 年 Rutland 和 Jameson[5]也由實驗觀察而得，尾隨液滴的大 小的確隨著波數的減小而增大。 1977 年 Pimbley 和 Lee[6]藉由實驗觀察尾隨液滴的形成，他們由 實驗得知，連接液滴的絲狀液流（Ligament）的分離處，可能發生在 前端、後端或前後端同時發生，絲狀液流可能是從下游處先分離，或 為從上游處分離；另外生成的尾隨液滴可能會與鄰近的主液滴合併， 如果絲狀液流同時在前後端分離，此時生成的尾隨液滴其速度與主液 滴相同，因此不會有合併的情形發生；另外當振幅小的時候，尾隨液 滴會與後方的液滴結合，而振幅大的時候則與前方液滴結合，而噴墨 列印技術所希望的現象，即是尾隨液滴不與主液滴結合的情形。

圖 2-2 尾隨液滴往前端合併

在噴墨的模擬上，歷史如下: 1974 年 Lee[7]首先提出了無黏度的模型，用以描述液體噴流形成 液 滴(而在多年後，考慮黏度效應的模型才被提出，如 1994 年 Eggers 和Dupont[8]) 。 1981 年 Kyser[9]等人提出半經驗式的模型(DOD 噴流的型式) 。 1982 年 Fromm[10]使用了軸對稱的Navier-Stokes 分析。

1984 年 Fromm[11]說明了當雷諾數（Reynolds Number）和韋伯數

（Weber Number）比值小時，操作條件必須使用較強的壓力。 1984 年 Bogy 和 Talke[12]利用實驗和模擬進行壓電式的研究，他 們改變不同的凹槽（Cavity）長度，比較液滴的生成時間、速度等。 1985 年 Allen 等人[13]利用軸對稱方式計算氣泡成長和液滴的生 成。 1992 年 Asai[14]利用三維的Navier-Stokes 方程式模擬氣泡的增長 和液滴射出的行為，並討論墨水的性質對液滴大小和速度的影響。 1998 年 Ho 等人[15.16]提出了新的設計概念來消除跟隨液滴，並增 加液滴的射出頻率。 1999 年 Chen[17]等人，由模擬與實驗研究壓電式噴墨頭，討論壓 力變化和液滴生成情形。

2002 年 Liou[18]等人結合VOF（Volume of Fluid）進行三維模擬，

研究液滴形成及其界面變化；2002 年 Brünahl[19]等人進行壓電式剪切 型驅動器（Actuator）的研究，討論不同條件對液滴的影響。 改善噴墨列印技術列印品質及速度，一直是我們追求的目標，因 此對於縮小液滴體積，增加噴墨頭上的噴嘴數目，增加射出速率，降 低尾隨液滴干擾問題等，皆有相當多的研究投入；然而當尺寸小至奈 米尺度時，傳統連續力學的適用性是一個相當大的問號，而且在如此

小的尺度下，整個噴流所表現行為，可能與連續力學下有所差異，因 此對於奈米噴流的研究，勢必改用其它方法進行。

利用分子模擬進行液體噴流的研究方面，分別如下:

1993 年 Koplik 與 Banavar[20]進行界面斷離（Rupture）的研究，

為他們研究液柱斷離形成液滴的過程圖，系統由一開始的液柱逐漸斷 離，然後慢慢調整，最後成為圓球；他們改變不同的液柱半徑，並將 液柱斷離的時間與線性理論做比較，得到了不錯的結果。

1998 年 Kawano[21]也進行了液柱在兩相斷離的研究，如圖 2.5， 他加大了系統邊長與液柱半徑的比值，所生成的液滴顆數並不像 Koplik 與 Banavar 僅有單顆，而下圖所生成的顆數較多，這是因為上 圖的液柱半徑較大；而Kawano 也將結果與 Rayleigh 的線性理論相比 較，且得到相當不錯的結果。 圖 2-5 液柱斷離過程圖（Kawano） 2000 年 Moseler 和 Landman[22]進行奈米噴流的研究，如圖2.6、 2.7，他們成功的模擬連續噴流現象，進而討論整個液流的形成、穩 定性及分離的現象，並且利用連續力學的概念，將 LE（Lubrication

Equation ） 加 上 一 修 正 項 項 成 為 SLE （ Stochastic Lubrication Equation），且由所得的 SLE 與分子模擬比較，得到與分子模擬相當

吻合的結果。Landman 等人所模擬的奈米噴流情形，這是考慮分子與

噴嘴有潤溼（Wetting）現象，由圖中可看出液柱及液滴皆出現非對稱 情形，且伴隨著蒸發的分子；而在傳統連續力學模擬中，一般常見的

方式是將Navier-Stokes 方程式簡化成一維軸對稱，因此在噴流過程中 有些現象無法忠實的呈現，雖然有人採用三維的模擬方式，可有效的 描述非對稱情況，但是當尺度小至奈米級時，對於分子間作用力所產 生的現象，仍然無法用巨觀的Navier-Stokes 方程式來描述，在液滴形 成前，在液流頸部（Neck）MD 是呈現雙錐形，利用連續力學的 LE 所得到的結果呈現細長的液流，SLE 模擬所得到的結果則與 MD 相當 穩合。 圖2-6 奈米噴流過程圖 圖2-7 液體噴流預測情形

2001 年 Goto[23]等人經由實驗與模擬進行研究，他們主要是利用 雷射熔射法（Laser Ablation），將一種有機分子注入另一種有機分子

表面，如何成功的將分子射出，如圖2.8、2.9。

圖2-8 雷射熔射法情況一

而在奈米噴流的實驗上，Voigt[24~27]等人做了一系列的研究，他 們由製作的奈米噴嘴，進行不同條件的奈米級蝕刻研究，他們發現， 當基板與噴嘴的距離小於噴嘴的半徑時，蝕刻的寬度會趨近於噴嘴的 直徑，亦即可產生約100 nm 的高解析度；而因為單一噴嘴蝕刻速度 太慢，所以為了克服這個問題，他們同時也發展了奈米噴嘴陣列，如 圖2.10。 圖 2-10 奈米噴嘴陣列 而奈米噴嘴製作的過程，利用聚焦式離子束（Focused Ion-Beam） 在金字塔型的矽後方鑽孔，而噴嘴的幾何主要受到離子束流動分佈 （Profile）的影響，因此經由適當的調整控制，可得到理想的噴嘴幾 何。 奈米噴流的應用相當的廣泛，而當小至奈米級時，利用噴流的特 性，可對MEMS 或 NEMS 的設備進行偵測和修補，或者製造精細的 奈米結構，例如EUV（Extreme Ultraviolet）光罩，如圖 2.11，另外也 可對表面進行平坦化，此有助微光學元件的製作，而對於生物細胞或 分子，也有其相當的應用，所以對於奈米噴流的研究，將有助於未來 在各領域的發展。

圖 2-11 奈米噴嘴的製造 2002 年 Tanner[30]模擬 C30的同分異構物在奈米尺度下剪切流場 (Shear Flow)中的行為，如圖 2.12。研究發現高分子之黏度與其支鏈 程度、第一正向應力有關，且在同一外加剪切應力下，支鏈越少的高 分子的滑動(Slip)程度越大。 圖2.12 剪切流場示意圖

2003年Xi-Jun Fan[31]利用分子動力研究奈米尺度下較複雜系統之 孔道流動，同時也利用有限元素法來做分析，如圖2.13。結果以分子 動力學模擬的結果會有渦流(Vortex)現象產生，並且在收縮管處有明 顯的滑動現象，但是以有限元素法模擬的結果則為無渦流的蠕流 (Creeping flow)。 圖2.13 收縮管流之模擬

2004年Te-Hua Fang[32]等人研究氬(Ar)在奈米噴流過程中孔洞及溫

度造成的影響，如圖2.14、2.15。根據分析，噴出的原子在溫度上升 時分布的更平均。當溫度下降，原子更輕易的集中在中間區域，當奈 米噴流孔洞更大時，噴出原子的量增加且更集中在中間區域，壓力分 布的現象可用分子動力學模擬得知。

圖 2.15 相同溫度不同管徑之噴流

另外分子動力學的發展史與電腦的發展史息息相關，Metropolis

在 1953 年第一個在 MANIAC 電腦上以蒙地卡羅法(MC)來模擬液體

。而後Alder 和 Wainwright[20]在1957 年提出硬球(Hard Sphere)碰撞

理論來模擬分子間的交互作用，為分子動力模擬(MD)的開端。 1964 年，Rahman[21]第一個應用真實的勢能模型來模擬液態氬(Ar) 。很快的在 1968 年 Berne[22]等人已經能模擬雙原子分子的行為，而 後1974 年 Rahman 和 Stillinger[23]完成水分子的模擬，如圖2.16。然 而當分子動力學架構完備，往後的進展更是快速。在 1977 年， McCammon[24]第一個完成模擬蛋白質的研究。之後分子動力學更廣泛 應用於各個領域之中。如蛋白質摺疊、DNA 結構、高分子物性、材 料性質…等等。

圖2.16 水分子與 DNA 的模擬

2002年，Masakazu Matsumoto、Shinji Saito 、 Iwao Ohmine[41]

應用真實的勢能模型來模擬液態水分子的結冰模擬，如圖2.17。

(1) (2) (3) (4) (5) (6) (7) (8) (1) (2) (3) (4) (5) (6) (7) (8) 2006 年清華大學林俊儀 [50]以分子動力學模擬奈米壓印製程 (Nano imprinting)。圖 2.25 顯示高分子在壓印加工時，最大壓力及最 大剪切應力出現在模具轉角處。 圖2.18 奈米壓印製程模擬[50] 2006 年交通大學曾煥錩[51][47]以分子動力學模擬奈米尺度下高分 子之擠出流場(Extrusion Flow)。圖 2.19 顯示擠出過程中有模口膨脹之 現象產生。 圖 2.19 奈米擠出流場模擬[47][51]

2006 年清華大學戴啟夫[49]以分子動力學模擬奈米尺度下固態之 聚乙烯高分子受到兩側牆面的剪切而產生的碎裂情形。圖 2.20 之虛 線為初始狀態的速度分佈，實線為經過一段時間後的速度分佈。 圖2.20 奈米固態薄膜受剪切力碎裂之模擬[49] 圖 2-21[52]為以分子動力學模擬高分子在剪切流場中的行為， 在得出分子鏈的基本相關數據後，代入性質方程式後得到相關連續力 學性質，例如：黏度、彈性係數、彈性模組等，並探討在剪切流道中 的現象。 圖 2-21 模擬收縮擴張流埸[52] Z/σ

2006 年的戴啟夫[50]

，以分子動力學模擬 PE 高分子的玻璃轉 移溫度，在 x、y、z 三個方向採用週期性邊界來模擬塊狀巨觀的情形，

圖 2-22[50]即為其模擬塊狀性質的示意圖，其測得數值符合實驗數據。

第三章、研究方法

3.1 分子動力學理論

分子動力學是一門藉由計算分子間受力情形，進而模擬出分子運 動軌跡的科學，可用來模擬真實狀態的分子行為，並且藉此1 來求得 系統中的各種特殊性質；整個分子間的作用力計算方式採用古典力學 (牛頓運動定律)的物理行為。針對分子動力學理論和研究方法分別敘 述如下。

3.1.1 分子動力學基本假設

分子動力學是以牛頓力學作為基本架構，理論計算中做了以下假 設： 一、所有分子皆遵守古典牛頓力學的運動定律 將分子間交互的作用力單純的簡化在勢能方程式上，忽略 其他效應的作用現象。 二、粒子的交互作用滿足疊加原理 意即將所有對分子作用的各種效應累加在一起，包括量子 效應和多體作用力，將其整合於一勢能函式上，因此，當以量 子力學或其他散射實驗測量方法所求得的分子勢能函式，即總 括了所有的分子受力效應，之後再利用古典力學來建立分子的 運動方程式。 三、假設每個分子均為球型，並以質心作計算

3.1.2 分子動力學模擬流程架構

根據上節所述的基本假設，再搭配合適的數值演算法與分子勢 能，就可以進行分子動力學模擬計算，而完整的分子動力學中有幾個 主要基本架構，可以分為五個部分： 一、初始化條件 分子動力學模擬的起始，需先給定模擬系統的初始狀態，由於 系統內的分子間受力是根據牛頓運動定律求得，因此在系統沒 有特殊需求的情況下，在初始設定只需給定初始位置r與初始 速度v。 二、系統控制 根據模擬的系統類型，需針對此模擬系統狀態條件分別作控 制，控制的部分包括溫度(T)控制、壓力(P)控制、體積(V)控制、 邊界條件設定等。 三、分子受力的計算 分子受力的計算為分子動力學的核心，根據牛頓運動定律並配 合適當的分子勢能函式，計算出系統分子在當時所受加速度 a，以便對下個時間分子運動軌跡作預測。 四、位置的預測 將已知的分子位置r和速度v，再加上計算所得的加速度a，配 合適當的數值演算法，即可求得新的步進時間(Time Step)下， 分子的新位置與速度。 五、計算性質參數 在計算上列分子運動行為時，可以得到每個步進時間下每個分 子的r、v、a等相關數值，將這些數值經過統計平均處理過後，

就可以將原本的微觀性質(Microscopic Properties)計算出巨觀性 質(Macroscopic Properties)如系統溫度、系統壓力、流體黏度、 應力張量、熱力學性質等。 整合以上每個基本架構，就可以構成完整分子動力學的運算過 程，其模擬流程圖如圖3.1 所示。 圖3.1 分子動力學模擬流程圖

3.1.3 分子動力學系統初始化

執行分子動力學模擬的第一個步驟，就是將系統內所包含的分子 進行初始化，主要處理在於設定分子的初始位置(r0)和初始速度(v0)， 模擬系統初始化設定— 模擬分子的位置r 初始化 模擬分子的速度v初始化 迴圈 依牛頓運動定律&分子勢能方程式— 模擬分子的所受加速度a 數值方法求解運動方程式— 求得模擬分子於(t +△t)時的位置r 求得模擬分子於(t +△t)時的速度v 結束判斷 根 據 模 擬 過 程 中 的 資 料 統 計 分 析 相 關 參數性質特性 未達最後時間 t = t +△t 是 否

由初始化的好壞將會影響到系統穩定與否，不佳初始化設定甚至會導 致模擬計算發散。

3.1.3.1 系統位置初始化

位置初始化通常會以真實物質最穩定的狀態來設定，如此一來， 不僅可以使系統較快達到穩定狀態，同時也可避免不必要的模擬錯 誤，因此，一般分子初始位置設定會依照物質的種類來區分，當固態 分子設定初始位置時通常會採用其個別晶體的結構，若為氣態或液態 時就會使用任意的晶體位置或亂數設定。 固態分子初始位置設定依照自然界的原子排向結晶狀況決定，在 自然界存在的結晶分別有三種堆積結構，其中分別所代表的結構有: 面心立方晶體結構(Face-Centered Cubic Crystal Structure) 、體心立方 晶體結構(Body-Centered Cubic Crystal Structure)、六方最密堆積晶體 結構(Hexagona Close-Packed Crystal Structure) 。

當所安排的分子呈現氣態或液態時，也可以使用上列的晶格排列 作為其初始位置，只需使系統進行適當的模擬時間讓粒子呈現較真實 的混亂分布，又或是初始就給定亂數的位置，可以根據Maxwell 分布 曲線來篩選亂數，但須注意兩粒子間的距離不能過近，否則會因受力 過大而產生外溢的現象。 此外，當模擬具有特殊鍵結的分子結構時，如:高分子鏈、蛋白質 結構等，通常會以平衡鍵長與平衡鍵角作為其起始結構，或是採用蒙 地卡羅(MC)方法製造出較穩定結構的初始位置。

3.1.3.2 系統速度初始化

執行分子動力學速度初始化時，最常使用的方法是採用亂數給定 速度，但為了使模擬初始合理化，通常會根據Maxwell 分布曲線來篩 選亂數的給定，並將亂數速度限定在一個速度範圍內，此外，如有特 殊的模擬系統需求時，將會另外考慮初始化動量歸零或是角動量歸 零。 1.初始動量歸零 先將所有分子的初始動量累加，再平均求得每個粒子的平均動 量，將每個分子減去此平均動量，則初始總動量即歸零。 i N k k k i i new i m m N m 1 ⎟/ ⎠ ⎞ ⎜ ⎝ ⎛ ₋ = v

∑

v v (3 – 1) 2.初始角動量歸零 類似於動量歸零的方法，將所有粒子的角動量累加後平均，再將 各個粒子的角動量減去此角動量則可使角動量歸零。

3.1.4 分子動力學系統控制

分子動力學模擬依據所模擬的系統，因應所採用的邊界條件、分 子控制方法不同，而有不同的處理方式，根據模擬系統種類可區分為 三類： 一、微正則系綜(Microcanonical Ensemble) 為 NVE 系統，模擬過程中，粒子數(N)、系統體積(V)以及系統總能 量(E)為一定值。

二、正則系綜(Canonical Ensemble)： 為NVT 系統，模擬過程中，粒子數(N)、系統體積(V)以及系統溫度(T) 為一定值。 三、NPT 系統： 模擬過程中，粒子數(N)、系統溫度(T)以及系統壓力(P)為一定值。 分別根據模擬系統需求作系統控制的動作，以下將針對如何設定系 統溫度控制、分子動力學求解運動方程式以及模擬系統時間間隔分別 作討論。

3.1.4.1 系統溫度控制

於恆溫的模擬系統中，需要使整體系統的溫度(T)保持一定值，由 於溫度乃是系統分子動能總合的表現，因此透過分子速度的調整，即 可維持系統的溫度；調整粒子速度的方法有下列幾種：

一、Simple Velocity Scaling

直接根據氣體動力學理論和動能方程式可得式(3 – 2)和式(3 – 3)： kT U_k 2 3 = (3 – 2)

∑

⋅ = i i i i k m N U v v 2 1 _{(3 – 3)} 其中Uk表示分子動能之總和，N 為系統總粒子數，m 表示分子的 質量，v 表示分子運動速度，k 為 Boltzmann’s 常數，結合二式則可以 得到：

∑

⋅ = i i i i m Nk T v v 3 1 _{(3 – 4)}

意即溫度即分子的速度平方總和，由式(3-4)就可以計算系統真正 的溫度 TA，若我們設定的溫度為 TD，則分子的速度就可由式(3-5)直 接進行調整： A D new T T v v = (3 – 5) 若為流動系統時，則分子動能總和需先減去流場之平均速度Ui， 才能代表系統分子之真實動能總和。

(

)

∑

− = i i i i k m N U 2 2 1 U v (3 – 6) 二、Nosé-Hoover Thermostat W.G. Hoover[30]於1986 年根據 S. Nosé[31]所提出的修正牛頓運動方 程式作進一步的改寫，因此，此方程式稱為Nosé-Hoover Thermostat： i i i m p v = (3 – 7) i i i F p p = −ζ ⋅ (3 – 8) ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ − =

∑

i D i i gkT m Q dt dζ 1 p2 _{(3 – 9)} 其中ζ 稱為摩擦係數(Friction Coefficient)，g 為系統自由度，Q為 溫控質量Q= gkT_Dτ2，τ 為控溫的特徵時間，利用此三式，可修正每個 分子的受力關係，進而達到修正系統的速度來控制系統的溫度， Nosé-Hoover Thermostat 是 NVT 系統中最常被使用的方法。 此方法的優點在於能夠維持系統能量恆定，較能描述真實分子運 動之行為。

3.1.4.2 模擬時間步進

時間步進(Time Step)為每次預測分子位置的時間差，時間步進選 取的大小將直接影響到模擬結果的準確性，其大小可以視系統變化程 度、分子間作用力的變化而有所增減，時間步進的選取也會隨著系統 密度增加、溫度上升、壓力增加而縮小。 使用時間步進過大會使模擬失真，時間步進過小則會拉長計算時 間，時間步進大小的決定通常由系統種類而定，一般可參考文獻所使 用之時間步進，或藉由嘗試修改時間步進來觀察其變化結果加以修 正。 在過去相關的分子動力學文獻中[28]，也曾經有提出使用多重時間

步進(Multiple Time Step)分別對不同作用力使用不同的時間步進，當 計算作用範圍大但作用力較小的作用力如凡得瓦爾力，可使用較長的 時間步進，而作用範圍小且作用力大的作用力如鍵的拉伸、彎曲、扭 轉，則採用較小的時間步進。

3.1.5 分子動力學求解運動方程式

分子動力學中的求解運動方程式，主要是利用泰勒展開式將每個 分子的位置和速度進行展開： ) ( ) ( 6 1 ) ( 2 1 ) ( ) ( ) ( 3 4 3 3 2 2 2 t O t dt t d t dt t d t dt t d t t t+∆ =r + r ∆ + r ∆ + r ∆ + ∆ r (3 – 10) ) ( ) ( 6 1 ) ( 2 1 ) ( ) ( ) ( 3 4 3 3 2 2 2 t O t dt t d t dt t d t dt t d t t t+∆ =v + v ∆ + v ∆ + v ∆ + ∆ v (3 – 11)

再根據此二式，推導出Gear Predictor-Corrector Algorithm 有限差分型式的演算方法： 由A. Rahman[29]首先將 Predictor-Corrector 使用在分子動力學上， Gear[30]再將此方法整理為一個完整方法；此求解運動方程式方法的處 理過程分成三個部分，分別為預測(Prediction)、計算(Evaluation)和修 正(Correction)。 首先利用已知的數值計算出下個時間間隔分子的位置與其位置 的微分值，第二個步驟求出預測位置的分子所受的加速度，最後再將 此加速度和上個時間間隔的加速度差值修正預測，求出正確的粒子位 置和速度，以下針對五階Predict-Corrector 方法作展開。 1.預測：首先使用泰勒展開式對位置、速度、加速度以及更多的微分 項作預測，其展開式如下： ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ ∆ ∆ ∆ ∆ ∆ ∆ ∆ ∆ ∆ ∆ ∆ ∆ ∆ ∆ ∆ = ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ ∆ + ∆ + ∆ + ∆ + ∆ + ∆ + ) ( ) ( ) ( ) ( ) ( ) ( 1 0 0 0 0 0 1 0 0 0 0 1 0 0 0 1 0 0 1 0 1 ) ( ) ( ) ( ) ( ) ( ) ( 2 ! 21 3 ! 31 2 ! 21 4 ! 41 3 ! 31 2 ! 21 5 ! 51 4 ! 41 3 ! 31 2 ! 21 t t t t t t t t t t t t t t t t t t t t t t t t t t t t t t t t t V IV III II V IV III II r r r r v r r r r r v r (3 – 12) 上標羅馬數字代表微分次數。 2.計算：由牛頓第二運動定律求出在新預測位置每個分子所受的力， 並轉換為分子加速度。 r F( ) ( ( ))ˆ r t t r U t t ∂ ∆ + ∂ − = ∆ + (3 – 13) I i m t t t t ) ( ) ( +∆ = F +∆ a (3 – 14) 3.修正：將步驟 2 計算出的∆rII(t+∆t)與步驟 3 中的a_i(t+∆t)相減求出

差值： )] ( ) ( [ ) (t t _i t t II t t II +∆ = +∆ − +∆ ∆r a r (3 – 15) II t r R≡ ∆ ∆ ∆ 2 2 (3 – 16) 將原先預測值做修正： R r r r r v r r r r r v r ∆ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ + ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ ∆ ∆ ∆ ∆ ∆ = ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ ∆ ∆ ∆ ∆ ∆ 5 4 3 2 1 0 5 ! 51 4 ! 41 3 ! 31 2 2 1 5 ! 51 4 ! 41 3 ! 31 2 2 1 α α α α α α t t t t t t t t t t V IV III II V IV III II (3 – 17) 於式(3 -28)中的α0、α1、α2、α3、α4、α5為Gear Predict-Corrector 的常數值，會隨著演算法計算階數(q)不同而改變參數值，其數值如表 3.1 所示。

表3.1 三～五階 Gear Predict-Corrector Algorithm 參數表

αi q = 3 q = 4 q = 5 α0 61 12019 163 α1 65 43 360251 α2 1 1 1 α3 31 21 1811 α4 － 121 61 α5 － － 601

3.1.6 分子動力學簡化方法

分子動力學的過程相當耗費計算時間，主要的原因在於模擬系統

的大小和模擬時間的長短，可以預期的是分子數目通常都需要龐大的 數量，而模擬的時間間隔數通常也在數十萬以上，因此適當的簡化計 算過程為分子動力學中重要的一環，否則很難將模擬時間控制在合理 的範圍內，如此一來，即使能正確模擬出重要的現象，但無法在合理 時間內完成模擬也是不合實際的。 由於過去的的電腦運算速率相較於現今明顯不足，因此，在早期 就有許多簡化計算的方法提供分子動力學，以下就常見的簡化方法作 個別介紹。

3.1.6.1 分子勢能截斷法(Cutoff Distance)

不論分子間使用的勢能方程式為何種類型，當分子與分子間的距 離 超 過 某 一 段 距 離 後 ， 分 子 間 的 作 用 力 將 會 變 的 極 小 ， 以 Lennard-Jones 勢能方程式為例，其勢能方程式對距離的關係圖(圖 3.2)： 因此可以設定一個截斷距離(Cutoff Distance，rc)，當分子間的距 離大於此距離時，直接假設此兩分子間的作用力為零，如此一來，就 可節省要計算遠距離的分子作用力，且因於截斷距離外之分子作用力 影響微乎其微，對分子預測幾乎不造成影響，卻由能有效簡化龐大的 計算量；大部分的分子動力學文獻，都會採用分子勢能截斷法來進行 模擬。 截斷距離長度的選擇，是由勢能函式的變化程度決定的，如模擬 系統僅包含簡單分子，且不考慮靜電效應與庫倫力時，通常所使用的 截斷距離約在 2.5σ~3.0σ 之間，σ 為兩分子間平衡距離；反之，若考 慮具有電荷的分子行為時，則截斷距離會選擇在 5σ 以上甚至超過

10σ。

distance between two atom, r*

0.5 1 1.5 2 2.5 3 -3 -2 -1 0 1 2 3 repulsive

force Attractive force

圖3.2 Lennard-Jones 勢能方程式及作用力

3.1.6.2Verlet 鄰近列表法(Neighbor List)

分子動力學模擬最耗時的部分在於判斷兩分子間距離，當系統特 別龐大時，此種判斷會更為耗時，即便是分子勢能截斷法，也需先求 得分子間距離，因此，運算量會隨著分子數目的平方關係成長，一旦 系統內分子數過多時，則會大量拖累運算速度；於是在1967 年 Verlet 提出鄰近列表法[9]來減少判斷分子間距離之運算。 鄰近列表法的基本觀念，在系統中的任一分子，其周圍的分子於 短時間內變動並不會太大，意即於此段時間內，此分子僅需和周圍鄰 近固定的分子計算作用力，可以不用對整個系統的分子作計算；當此 方法使用在越大的系統時，加速的效果更為顯著。 Potential U*(r) Force F*(r)

鄰近列表法的使用方法，先將分子i 列表半徑(List Distance ,rl)內 的分子都列入清單中，在接下來的時間間隔中，分子i 只需對清單內 的分子作距離計算，持續到下一次的更新列表，而鄰近列表清單只需 要數個或數十個時間間隔後再予與更新。 其中列表半徑和清單列表的頻率是相互影響的，選擇的決定依據 在於，位於列表半徑外的分子能否在下次清單更新前，進入影響分子 的截斷距離內，一般而言，當分子的移動速度越快，相對的列表更新 頻率增高或者是列表半徑要增加。 圖3.3 為鄰近列表法的簡易的示意圖： 圖3.3 鄰近列表法示意圖 圖3.3 顯示以 1 號分子為中心，2 到 13 號為其周圍的分子，虛線 圓圈為1 號分子截斷半徑的範圍，實線圓圈則為列表半徑的範圍，每 個時間間隔下，1 號分子要計算分子間距離時僅需考慮到 2 到 9 號分 子，因為10 和 13 號分子在下次列表前，都無法進入截斷距離的範圍 內，故不需要列入計算 8 7 3 1 2 6 10 11 4 5 9 13 12

3.1.6.3 Cell-linking 列表法

鄰近列表法可以有效的減少每個分子對距離較遠分子作用力的 計算，但是每當進行列表更新時，仍舊是個耗時的龐大計算，因此隨 後Verlet 又提出了 Cell-linking 列表法，此方法使運算量僅與分子數目 成一次方關係，大量減少運算量，並且可與鄰近列表法做結合，更有 效率的節省計算時所花費的時間。 此 方 法 的 概 念 是 先 將 模 擬 的 系 統 切 割 成 若 干 等 大 小 的 方 格 (Cell)，如圖 3.4 所示。 圖3.4 Cell-linking 列表法示意圖 判斷出系統中每個分子所處的方格，當欲計算系統某一個分子受 力情形時，只需找尋在分子所屬方格周圍附近，未超過參考距離的數 個方格內的分子來判斷即可，不需要對整個系統的分子個別做搜尋判 斷，而考慮使用的距離應為截斷距離，但絕大部分此方法都會和鄰近 列表法同時使用，則參考距離就採用列表半徑。 在圖 3.3 中，以一個 2D 系統為例，當要計算位於深色方格中的 分子受力時，只需要考慮以深色方格為中心，且未超過考慮距離的九

個方格中的分子，同理當使用在 3D 系統時，每一個方格中的分子也 只需要考慮周圍 27 個方格中的分子，相較原先需要考慮整個系統其 他分子的方法簡化了許多。 此外，切割方格的大小也會影響其運算效率，當切割的越大時多 餘的判斷會增加，而相反的切割的過小，也會造成找尋的方格數增 加，因此選定一個適當的方格長度也是相當重要的，一般而言，選定 列表距離的一半通常會有較佳的效率。

3.1.6.4 列表式勢能(Tabulated Potential)

計算分子受力為分子動力學的主要核心，計算的方法是將兩分子 間的距離帶入微分後的勢能方程式得到受力情形，而分子間距離通常 也只在零到截斷距離的長度而已，出現相同分子距離的機會相當大， 且受力情形和分子間距離的變化並不是太大，因此在 1999 年由 D. Wolff 等人提出以列表式勢能[32]來簡化計算的步驟。 在開始分子模擬前，先將分子間距與此間距下的作用力計算出並 存入記憶體中，而後只需將分子間距離算出後直接找出對應的作用 力，不需經過函式的計算來節省時間；因此，D. Wolff 等人提出利用 此方法將勢能方程式切割成好幾個等分，製作成階梯供查詢，圖 3.5 為Morse 勢能方程式分割為 50 個等分後製作而成的列表式勢能圖。

r (Angstroms) Mo rs e E n e rg y 0 1 2 3 4 5 6 7 8 -100 -50 0 50 100 150 200 圖3.5 Morse 列表式勢能圖 使用此列表式勢能需注意於作用力分割的精細度，當分割階梯越 多模擬會越準確，但是使用的記憶體相對會增加，容易造成電腦的負 擔；分割的階梯越少其階層的變化越大，容易造成誤差，因此階梯分 割的精細度需根據勢能方程式的變化以及系統平均分子間距離的情 形來做調整。 列表式勢能方法通常只用在相當複雜的勢能函式，如多體勢能或 Morse 等勢能方程式，在計算過程中會使用到需大量浮點運算的算式 如指數運算，其求解時間明顯會比列表後直接取值慢，使用此方法可 以 節 省 這 類 運 算 的 時 間 ， 但 使 用 在 簡 單 運 算 的 勢 能 方 程 式 ， 如 Lennard-Jones 勢能方程式，則沒有太明顯的加速效果。

3.1.7 分子動力學性質計算

分子動力學可以準確利用牛頓定律模擬分子行為狀態，但所得資 訊原先僅有表觀狀態，即分子運動軌跡狀態而非真實世界中所觀測到 的性質，因此需要經過統計力學相關推導才能引入巨觀的性質行為， 所包含的相關函式如下： 一、基本性質 1.溫度(Temperature) 在微觀角度下，我們所得到的資訊並非溫度，而是由分子運動所 攜帶的能量，因此如果要得知模擬系統的溫度，就需要藉由計算所得 的分子速度分析而得： > < = E_k N T 3 2 _{(3 – 18)} 2 2 1 mv E_k = (3 – 19) 2.能量(Potential energy) 由已知的分子間勢能可推導出系統能量之變化，由於分子勢能計 算都會採用截斷距離來處理，因此還需要增加一項長距離(long range) 勢能進行修正，式(3 – 22)為相同分子於截斷半徑外的能量總和： LR cutoff total U U U = + (3 – 20)

∑∑

> = N i N i j ij cutoff U r U ( ) (3 – 21) 3 6 3 ₃ 8 1 3 1 3 8 c c c LR r r r U πρ _⎟⎟≈− πρ ⎠ ⎞ ⎜⎜ ⎝ ⎛ − = (3 – 22) 3.壓力(Pressure) 分子尺度所造成的壓力可分為兩個部分，第一部份是分子通過一

平面的動量通量(Momentum Flux)，第二個部分則是分子間引力造成 的能量場差，因此壓力表示式可寫為：

∑

∑∑

> ⋅ + ⋅ − = N i N i N i j ij ij i i i m V P v v r F 3 1 (3 – 23) 4.數目密度分布(Local Number Density)

用來形容系統分子間的分佈距離關係，可用來作為判斷系統分子 相態的工具： ) ∆ , ( ∆ ) ( ) , ( ) , ( ) ( r r V r r r dr r V dr r N r N i i

∑

− = = δ ρ (3 – 24) N 為以分子為中心，r為半徑厚度，∆r的圓球薄殼中的分子數目， V 為此薄殼圓球體積。

3-2 勢能函數選擇

3.2.1 分子間勢能:

對於分子間的勢能函數，在本研究中採用Lennard-Jones 勢能

( )

⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ − ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ = 6 12 4 ij ij ij LJ r r r U

ε

σ

σ

_{(3 – 25)} 式中σ為兩分子間勢能等於零的距離，ε為兩分子間的最低能量，上 述兩參數皆隨分子種類不同而不同，而rij為兩分子間的距離。為了節 省計算量，因此在模擬時使用截斷半徑，所以勢能改寫為：

( )

⎪ ⎪ ⎩ ⎪⎪ ⎨ ⎧ > ≤ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ − ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ = c c ij ij ij LJ r r r r r r r U 0 4 6 12

σ

σ

ε

(3 – 26) 因為採用截斷半徑，所以勢能函數並不是呈現平滑曲線，這可能會導 致總能量無法守恆，為了避免這種情況發生，因此對整個函數做平 移，而上式改寫為：

( )

( )

( ) (

)

⎩ ⎨ ⎧ > ≤ ∆ − − − = c c c c LJ ij LJ ij Shift r r r r F r r r U r U r U 0 (3 – 27) ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ − ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ − = ∆ 7 13 2 24 c c r r F

σ

σ

σ

ε

(3 – 28) 式中rc為截斷半徑。 對於金屬原子的模擬，一般都採用多體勢能，由相關的文獻可 知，在各種情況下都有準確的描述，但是因為多體勢能在計算上相當 耗費時間，且就本研究的系統中，金屬原子皆處於晶格位置掁盪，因 此我們採用較簡單的Morse 勢能來描述，其勢能函數如下：

( )

[

2 (r r0)

2

(r r0)

]

ij Morse ij ij

e

e

D

r

U

=

− α −

−

−α − _{(3 – 29)} 在計算分子間作用力時，系統中可能包含不同的原子或分子，因

此計算時，分子交互間的參數由Lorentz-Berthelot Mixing Rule 求得，

(

a b

)

ab

σ

σ

σ

=

+

2

1

(3 – 30) b a ab

ε

ε

ε

=

(3 – 31) 本研究所使用的 Lennard-Jones 參數如下表所示，其中氫原子與氧原

子間的能量是由Materials Studio 而得，而 Lennard-Jones 勢能採用的

截斷半徑為3.0σCH2，Morse 勢能的截斷半徑為 2.7σCH2。 non-bonding potential 參數如下： 非鍵結的分子間交互影響以L-J 分子勢能模型描述： ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ − ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ = 6 12 4 ij ij ij ij ij r r E ε σ σ (3 – 32)

3.2.2 分子內勢能數據:

本研究所模擬的系統中，因為包含氫分子，氧分子以及 CH2 bead，所以需要考慮分子內勢能，其中包含拉伸鍵能（Bond-Stretching

Energy）

、

彎曲鍵能（Bond-Bending Energy）及扭轉鍵能（Bond-Torsion

Energy）等，函式分別如下：

1. bond stretching potential：

律來描述其行為。其中式中ks 為伸縮常數，req 為平衡鍵長。拉伸勢 能以類似彈簧的簡諧運動描述： 2 0) (r r k U_s = _s − (3 – 33)

2.bond bending potential：

鍵角的彎曲位勢能亦相似虎克定律，不同於是描述鍵結間鍵角 的 簡協運動，其中式中 ks 為鍵角彎曲常數，θ 為三個鄰近原子 間的鍵角，θ0為平衡鍵角。兩相鄰鍵結會有穩定鍵角產生，使用類 虎克勢能描述： 2 0) (θ−θ = b b k U (3–34) 3. Electrostatic energy ij j i j i ele r e q q U 2

∑

∑

= (3–35)

第四章 模擬系統

4.1 純水物性的模擬系統

一．模擬系統介紹 純水物性的測量有許多的方法，其中的 PVT 性質是最常見 的方式，所以在此我們就使用分子動力學去模擬分子的 PVT 性質,與 文獻相互驗證。 我們選擇與文獻相同條件的模擬系統來模擬水分子,以便數 據比較,系統如下,在一正立方體的系統中，保持固定的水分子個數 (216 個水分子)，所以系統的總粒子數不改變，僅改變不同組別的系 統體積，以得到不同的密度，再改變系統的溫度，即可得到不同密度 的壓力作圖。 二．模擬方法 首先我們為了以少量的粒子去模擬巨觀的性質，我們在三個 方向(X,Y,Z)都使用了週期性邊界，以達到在各個方向都是無限大的情 形，以得到類似巨觀的性質。 系統為 PVT 系統，所以系統中的粒子總數，一直保持恒定。 系統控溫則是採用 Nose-Hoovêr Thermostate，使系統的溫度維持不 變。在預測粒子運動軌跡的方程式則是用五階Gear Predictor-Corrector 方程式。