以分子動力學研究高分子溶液之結構性質與奈米噴流行為

國

立

交

通

大

學

應用化學系

碩

士

論

文

以分子動力學研究高分子溶液之結構性質與奈米噴

流行為

Studies on the behavior of nanojet and structure properties of

polymer solution by molecular dynamics simulation

研 究 生：洪健超

指導教授：吳建興

Studies on the behavior of nanojet flow and structure properties of polymer solution by molecular dynamics simulation

研 究 生：洪健超 Student：Chien-Chao Hung 指導教授：吳建興 Advisor：Jiann-Shing Wu

指導教授：張榮語 Advisor：Rong-Yu Chang

國 立 交 通 大 學 應 用 化 學 系

碩 士 論 文

A Thesis

Submitted to Department of Applied Chemistry College of Science

National Chiao Tung University For the Degree of Master

in

Applied Chemistry July 2008

Hsinchu, Taiwan, Republic of China

以分子動力學研究高分子溶液之結構性質與奈

米噴流行為

學生：洪健超 指導教授： 吳建興 張榮語 國立交通大學應用化學系 摘 要 軟性電子具有可塑性、重量輕、耐衝擊等優點，潛藏龐大的商機。製 程上透過噴印技術將導電性材料「印刷」於基材上，因此微尺度的噴墨列 印製程技術受到各界的重視。此外印刷製程也可以運用於電子元件構裝 上，倘若能達到奈米尺度的噴印，電子元件的發展也將進入一個新的次元。 本文藉由分子動力學模擬研究高分子溶液奈米噴流技術，透過不同的操作 條件觀察噴流過程，選用的高分子噴流材料為聚氧乙烯溶液，採用 EA model 進行模擬，透過結構簡單、應用甚廣的水溶性聚氧乙烯幫助瞭解高分子溶 液的噴流過程。金屬勢能則選擇 Morse potential 來使用，水分子則使用 Simple point-charge (SPC/Fw) model 進行模擬，相較於其他水分子模型，由 於在性質方面的再現性與實驗結果相當一致，同時是水分子模型中最簡單 的三質心點類型，擁有較佳的計算效率，適合運用在大尺度的加工模擬的 研究。 研究的過程分為二個階段：建立高分子溶液，討論溶液結構並驗證模 擬之程式的正確性；隨後選用合適的金原子建立噴嘴，建立噴流裝置模型； 本研究著重於驅動力的影響，研究不同的推擠溶液的速度對噴流之行為造

Investigate behavior of nanojet and structure

properties of polymer solution by molecular

dynamics simulation

student：Chien – Chao Hung Advisors： Jiann-Shing Wu

Rong-Yu Chang

Institute of Applied Chemistry National Chiao Tung University

Abstract

Duo to plasticity, light weight, impact resistance and other advantages, flexible electronics have the huge potential business opportunities .Micro scale jet printing technology will play an important role since the conductive material can be printed on a substrate via this printing process .In addition , the printing process can also be applied to the packaging of electronic components. If the jet printing process can achieve nanometer scale, the development of electronic components will also enter a new territory .In this study a nano-jet printing technology with a polymer solution as a printing material was investigated using molecular dynamics and then observed the jet printing process by varying the operating conditions. Due to its simple structure and wide usage, the polyethylene oxide was utilized as a printing material to simulate and to understand the jet printing process, using explicit atoms model (EA) .Morse potential was chose to represent a potential model for the metal with regards to the water molecule. Simple point-charge (SPC/Fw) model was used to be a potential model, due to the reproducibility and consistency in the experimental

results for the properties of water. Meanwhile, this SPC/Fw model is the simplest three-sites model which possesses better computational efficiency, consequently leading to be suitably applied to the investigation for the large scale processing simulation.

The process of this study was divided into two steps. The first step was establish the polymer solution and then discussed the structure solution and identifies the correctness of the simulation program. Second step was to choose a proper metal atom to construct a nozzle and to establish a model for the system of jet flow. Finally, this study aimed to investigate the effect of driving force and to study the effect of plunger speed on the behavior of jet flow.

誌謝

兩年的碩士生活讓學生學習到許多東西也成長不少，當然也少不了長 期的日夜巔倒的生活以及連日的敖夜，讓整個身體變的很差。 在學習的過程中，感謝吳建興老師的指導，在教學的過程中十分用心、 內容詳盡 (至少比某些人來的好)。同時受教於清華大學張榮語教授，賴於 老師的教導與指引讓論文編寫的過程更加的順利、完成度更佳，由衷感謝。 感謝學長、姊們在程式及相關理論的指導，雖然我覺得八成好像都是 在自學。同學煜智、長福是這兩年中最重要的伙伴，不論是休閒育樂、學 業或研究的討論等，總是一起努力渡過。尤其是長福同學，因為實在太神 了，總是受到他的幫忙與指點，實在慶興能做為他的同學，不過也因為同 期有高手的存在讓自己壓力很大。另外，也感謝各位學弟們，總是提供趣 事，讓研究的過程不再沉悶艱辛。 同時也要感謝在交通大學的兩位同學及學弟們，多虧小童童經常陪我 看日出，老人沒事請吃飯，讓兩年的學習生活多采多姿。 最後，也要感謝我的父母及家人對我的支持與幫助，希望他們永保健 康，也祝福所有人未來一路順利。

目錄

摘要 ... i Abstract ... ii 誌謝 ... iv 目錄 ... v 圖目錄 ...viii 表目錄 ... xi 符號說明 ...xii 1. 緒論... 1 2. 研究背景與目的... 2 2.1. 研究背景 ... 2 2.2. 軟性電子與噴墨印刷製程 ... 3 2.2.1. 軟性電子... 3 2.3. 噴墨技術 ... 6 2.4. 研究動機與目的 ... 9 3. 文獻回顧... 11 3.1. 水分子勢能模型 ... 11 3.2. 高分子與高分子溶液 ... 16 3.3. 奈米噴流 ... 18 4. 研究方法... 24 4.1. 分子動力學理論 ... 24 4.2. 勢能模型 ... 24 4.2.1. 水分子勢能模型... 25 4.2.2. 聚氧乙烯勢能... 26 4.2.3. 金屬勢能... 27

4.3. 初始化位置及速度 ... 32 4.4. 系統類型與控制 ... 33 4.4.1. 系統體積控制... 33 4.4.2. 系統控溫... 34 4.5. 數值演算法 ... 36 4.6. 加速計算方法 ... 38 4.6.1. 牛頓第三運動定律... 38 4.6.2. 分子勢能截斷法... 38 4.6.3. 鄰近列表法... 41 4.6.4. Cell-linking method ... 42

4.6.5. 列表勢能法 (tabulated potential method) ... 43

4.7. long-range interaction... 44 4.7.1. Ewald sum... 45 4.8. 性質計算 ... 48 4.8.1. 溫度... 48 4.8.2. 壓力... 48 4.8.3. 徑向分佈函數... 49 5. 模擬系統... 50 5.1. 水與聚氧乙烯水溶液塊材系統 ... 50 5.2. 高分子溶液噴流系統 ... 54 6. 結果與討論... 57 6.1. 塊材水溶液結構性質 ... 57 6.1.1. 平衡鍵長與分子結構... 57 6.1.2. 平衡鍵角與分子結構... 58 6.1.3. SPC/Fw 勢能模型之分子結構 ... 58 6.2. 高分子溶液結構性質 ... 64 6.2.1. 溶液中 1, 2-二甲氧基乙烷與水分子結構 ... 64

6.2.2. 溶液中水結構... 65 6.3. 高分子溶液噴流特性 ... 68 6.3.1. 噴流模擬擠出現象... 69 6.3.2. 液滴自然斷裂與飛行現象 ... 89 6.3.3. 噴流過程中高分子之結構性質 ... 93 7. 結論與未來展望... 96 參考文獻 ... 98

圖目錄

圖 2-1 軟性電子紙的優勢分析5... 4 圖 2-2 富士通發展之軟性電子看 板6 ... 5 圖 2-3 富士通發展之電子閱讀 書，展現其輕量化的特色6... 5 圖 2-4 富士通軟性電子商品型 錄6 ... 5 圖 2-5 噴墨印表機技術之分類圖7... 7 圖 2-6 熱泡式頂噴型噴墨技術示意圖7 ... 8 圖 2-7 壓電式推擠型噴墨示意圖7... 8 圖 3-1 液柱斷離過程圖49 ... 18 圖 3-2 液柱斷離過程圖（Kawano）50... 19 圖 3-3 奈米噴流過程圖10 ... 20 圖 3-4 液體噴流預測情形10... 20 圖 3-5 雷射熔射法51 ... 21 圖 3-6 奈米噴嘴陣列53... 21 圖 3-7 Argon 於基材之噴流塗佈情形及傾斜噴嘴之塗佈情形 ... 22 圖 3-8 不同形狀之噴嘴示意圖15... 23

圖 4-1 勢能 Lennard-Jones 12-6 shifted potential 示意圖 ... 39

圖 4-2 使用 Switch function 之 Spohr 勢能示意圖57... 40

圖 4-3 鄰近列表法示意圖 ... 41

圖 4-4 Cell-Linking 法示意圖... 42

圖 4-5 Complementary error function 列表勢能示意圖 ... 43

圖 4-6 比較庫侖力與 Lennard-Jones 勢能降低至零的程度... 45

圖 4-7 系統鏡像延伸後所建立的球形系統 ... 46

圖 5-1 216 個水分子(H2O)之平衡態熱力學系統三維示意圖 ... 51

圖 5-2 45 條 1, 2-二甲氧基乙烷 (DME) 及 200 顆水分子 (XDME =0.18) 之平 衡態勢力學系統三維示意... 51

圖 5-3 塊材性質模擬流程圖 ... 53 圖 5-4 嘴口收縮角度為 135°，以金原子為材料所建立的噴嘴示意圖 ... 54 圖 5-5 嘴模型幾何示意圖 ... 55 圖 5-6 奈米噴流模擬流程圖 ... 56 圖 6-1 水分子 256 顆於不同平衡鍵長設定，氧-氧原子徑向分佈函數... 59 圖 6-2 文獻 18 中水分子與不同平衡鍵長之氧-氧原子徑向分佈函數... 59 圖 6-3 水分子 216 顆於不同平衡鍵長設定，氧-氫原子徑向分佈函數... 60 圖 6-4 文獻 18 中水分子與不同平衡鍵長之氧-氫原子徑向分佈函數... 60 圖 6-5 水分子 216 顆於不同平衡鍵角設定，氧-氧原子徑向分佈函數... 61 圖 6-6 文獻 18 中水分子與不同平衡鍵角之氧-氧原子徑向分佈函數... 61 圖 6-7 水分子 216 顆 SPC/Fw model 徑向分佈函數... 62 圖 6-8 文獻 18 中比較實驗結果與 SPC/Fw model 的徑向分佈函數... 63

圖 6-9 DME 水溶液 ODME-OWater及 ODME-HWater的徑向分佈函數 ... 66

圖 6-10 文獻 46 中 DME 水溶液的徑向分佈函數 ... 66

圖 6-11 DME 水溶液於 318 K。分別為 DME 莫耳分率(XDME)0.18 及 0.72 的 局部(local)與塊狀(bulk)水溶液密度比值 ... 67

圖 6-12 水獻 46 中，DME 莫耳分率(XDME)為 0.18 及 0.72 的局部(local)與塊 狀(bulk)水溶液密度比值 ... 67 圖 6-13 不同推板位置示意圖 ... 68 圖 6-14 推擠速度 40 m/s 聚氧乙烯水溶液噴流過程圖 ... 71 圖 6-15 推擠速度 80 m/s 聚氧乙烯水溶液噴流過程圖 ... 72 圖 6-16 推擠速度 120 m/s 聚氧乙烯水溶液噴流過程圖 ... 73 圖 6-17 推擠速度 160 m/s 聚氧乙烯水溶液噴流過程圖 ... 74 圖 6-18 推擠速度 200 m/s 聚氧乙烯水溶液噴流過程圖 ... 75 圖 6-19 不同推擠速度於 6/8 推板位置之二維密度分佈圖 ... 77 圖 6-20 不同推擠速度於 5/8 之二維密度分佈圖呈現分層現象 ... 78

圖 6-22 噴嘴區域幾何示意圖 ... 80 圖 6-23 不同推擠速度在 5/8 推板位置之密度與軸向座標關係圖 ... 80 圖 6-24 不同推擠速度在 5/8 推板位置之壓力與軸向座標關係圖 ... 81 圖 6-25 推擠速度 40 m/s 在不同推板位置之密度與 Z 軸關係圖... 82 圖 6-26 推擠速度 80 m/s 在不同推板位置之密度與 Z 軸關係圖... 82 圖 6-27 推擠速度 120 m/s 在不同推板位置之密度與 Z 軸關係圖... 83 圖 6-28 推擠速度 40 m/s 在不同推板位置之壓力與 Z 軸關係圖... 83 圖 6-29 推擠速度 80 m/s 在不同推板位置之壓力與 Z 軸關係圖... 84 圖 6-30 推擠速度 120 m/s 在不同推板位置之壓力與 Z 軸關係圖... 84 圖 6-31 不同推擠速度之質量通量與時間關係圖 ... 85 圖 6-32 不同推擠速度在 5/8 推板位置之溫度與軸向座標關係圖 ... 86 圖 6-33 推擠速度 40 m/s 在不同推板位置之溫度與 Z 軸關係圖... 86 圖 6-34 推擠速度 80 m/s 在不同推板位置之溫度與 Z 軸關係圖... 87 圖 6-35 推擠速度 120 m/s 在不同推板位置之溫度與 Z 軸關係圖... 87 圖 6-36 不同推擠速度於 5/8 之二維溫度分佈圖呈現分層現象 ... 88 圖 6-37 推擠速度 40 m/s 於溶液用盡後噴流自然飛行過程圖 ... 90 圖 6-38 推擠速度 80 m/s 於溶液用盡後噴流自然飛行過程圖 ... 91 圖 6-39 推擠速度 120 m/s 於溶液用盡後噴流自然飛行過程圖 ... 92 圖 6-40 推擠速度 80 m/s 不同區域與迴旋半徑及末端-末端距離關係圖 ... 93 圖 6-41 不同推擠速度於不同區域與迴旋半徑關係圖 ... 95 圖 6-42 不同推擠速度於不同區域與末端-末端距離關係圖... 95

表目錄

表 3.1 數種水分子勢能計算塊狀水溶液之性質18... 15

表 4.1 SPC/Fw model 勢能參數表18 ... 29

表 4.2 Poly (ethylene oxide) 與 1,2-dimethoxyethane 勢能參數表42, 43... 30

表 4.3 Au 相關勢能及參數56 ... 31

表 4.4 Polyethylene – Water 勢能參數表46... 31

表 4.5 三~五階 Gear Predict-Corrector Algorithm 參數 ... 37

表 5.1 研究塊材水溶液性質之系統操作條件... 52

符號說明

bond U ：分子鏈拉伸勢能 bend U ：分子鏈彎曲勢能 OH r ：分子鍵長 0 OH r _{：分子平衡鍵長} HOH ∠ θ _{：分子鍵角} 0 HOH ∠ θ _{：分子平衡鍵角} er Uint _{：分子間勢能} ra Uint _{：分子內勢能} σ ε、 _{：Lennard-Jones 12-6 勢能參數} i

q _{：i 粒子部分電符 (partial charge)}

k a _{：雙面角勢能參數} φ_{：雙面角角度} ij ij ij B C A、 、 ： Lennard-Jones exp-6 勢 能參數 D：Morse potential 勢能參數 ( )r S _{：開關函數(Switch function) 勢} 能修正式 on r _{：開關函數開啟半徑} off r ：開關函數關閉半徑 i m _{：i 粒子質量} i v _{：i 粒子速度} N_{：系統總粒子數} V _{：系統體積} T ：系統溫度 E：系統內能 L：系統邊長 a T _{：控溫系統真實溫度} d T _{：控溫系統設定溫度} ζ _{：摩擦係數} g_{：系統自由度} Q_{：溫控質量} τ _{：特徵時間} s U _{：平移勢能(shifted potential)} c r _{：截斷半徑} o ε _{：介電常數} s ε _{：相對介電常數} n：實部空間向量 k：虛部空間向量 B k _{：波茲曼常數} c U _{：long-range 修正勢能} ρ_：密度 k E _：動能 M ：統計之時間步進(time-steps)

1. 緒論

著眼於奈米的特殊性質，因此奈米科技被廣泛 (nanotechnology)的研究 著，其所涵蓋的領域甚廣，從基礎科學橫跨至應用科學，包括物理、化學、 材料、光電、生物及醫藥等1。由於奈米科技的蓬勃發展，因而發展許多相 關應用、技術並且產品化，所持有的價值與競爭力，對未來的發展一定是 不可限量，因此奈米加工製成技術成為十分重要的課題。分子動力學模擬 (Molecular dynamics simulation) 提供了強大的功能，在奈米尺度之下，透過 電腦模擬將能夠預測、解決許多未知的問題。 隨著科技的進步，奈米技術的知識已經建立了許多的應用技術，例如 奈米結晶材料、奈米粉體、複合材料、塗佈及奈米纜線的相關技術，運用 這些技術研發出許多應用產品。然而，不可忽視的是製程的重要性，在製 造運用奈米技術的商品時，如何選擇最佳的加工方式並且深入瞭解使其最 佳化，以達到快速量產、高良率及低成本等，是相當值得討論的方向。

奈米製程的模擬技術，其中奈米噴流加工製程 (Nanojet printing process) 或奈米噴墨列印製程 (Nano ink-jet printing process) 極具潛力，擁有許多的 優點且應用範圍相當廣泛。因此，預計透過分子動力學模擬深入探討奈米 噴流加工技術。

2. 研究背景與目的

2.1. 研究背景

目前全球工業發展趨向微縮化製程技術，奈米加工製程將取代現今的 加工技術，以提高生效能與更優秀的奈米元件。然而，對於奈米尺度之下 的知識與技術，仍然有許多不了解的問題及尚未發覺的現象，若要實際的 產品化並且大量的生產製造，必然有許多的問題需要克服。 在微奈米尺度之下，電腦模擬成為十分有利的工具。其中分子動力學 模擬在奈米材料及高分子材料之物性、化性及加工製程正快速的發展，其 具有以下優點： 1. 實際的實驗中往往需要價格十分昂貴的器材，電腦模擬具有相當低的成 本。 2. 奈米尺度的實驗不易測量，透過模擬的方式不但不受限於儀器、設備， 甚至能突破實驗無法做到的操作件，彌補實驗不足的地方。 3. 許多研究資料顯示，分子動力學模擬之結果與實驗觀察現象都能高度吻 合2、3_{，即便不能獲得精確值，在定性方面同樣能獲得良好的參考。} 4. 實驗的方式往往需要大量的時間、材料消耗以及人力，對於產品的開 發、應用將會是個障礙，甚至可能延誤商機。 近年來電腦算速度與記憶體容量大家提昇，隨著電腦效能的提昇，分 子動力學模擬將能夠越來越接近真實的奈米元件製程系統，甚至能夠解決 任何的問題。 奈米材料所具有的性質與巨觀尺度下的特性截然不同，因而巨觀尺度

透過連續力學所建立的物理現象並不適用於該範圍，必需以微觀的角度包 含分子、量子理論進行研究分析，才能掌握奈米科技的特性與實際的應用、 製程。分子動力學便扮演著重要的角色，該技術已經廣泛的運用在多種領 域的研究，如奈米高分子之流變學、奈米加工製程及奈米生物等領域。分 子動力學能幫助奈米尺度下的研究，不但能估算奈米材料的性質甚至是掌 握所有原子的動向，在製程方面也是十分有用的利器，提供製程的最佳化 設計。

2.2. 軟性電子與噴墨印刷製程

在眾多奈米加工製程技術之中，奈米奈米噴墨列印製程技術為高潛力 的奈米光電元件加工製程，尤其是以高分子為材料之光電元件，如奈米軟 性電子薄膜、奈米光子材料、奈米金屬導線、奈米矽基材等等方面的奈米 光電元件。 2.2.1. 軟性電子 光電元件中，尤其以軟性電子 (Flexible Electronics) 十分具有潛力。目 前全球電子研發機構都在從事革命性的「軟性電子」研究，工研院電子所 指出，2009 年全球軟性電子產值達 58 億美元，2012 年上看 235 億美元4_。 軟性電子技術可讓電子產品的應用領域進一步擴展，使電子產品具備輕薄 短小、符合人因工程、便利穿戴、柔軟安全、可自由捲曲的特性，發展出 適合行動的數位化多功能電子生活產品。 所謂的軟性電子從字面上解釋，即為對比於過去無法彎曲、型態固定 的電子裝置，由於採用有機材料以及薄塑膠或金屬薄片為基材，讓電子元

態的特性，在螢幕畫面不變時，是不耗電的，只有在更換螢幕時，才需要 耗電，與目前平面顯示器需要一直耗電的狀況是不相同的，也正因為具有 此特性，使得軟性顯示器的耗電量非常低，進而也減少了電池的更換頻率， 且觀看角度不受限制，更重要的是可以透過噴墨印刷製程將導電材料佈設 在可撓性基材上，採用低成本的卷對卷 (Roll to Roll) 方式來製作。 圖 2-1 軟性電子紙的優勢分析5

圖 2-2 富士通發展之軟性電子看

圖 2-3 富士通發展之電子閱讀

書，展現其輕量化的特色6

2.3. 噴墨技術

噴印加工製程的優點在於其為溶液加工製程(Solution process)、能精準 控制液滴的量與位置、生產速度快速、成本較其他製程低廉、列印圖案可 具有高解度、減少半導體微影製程的次數、不接觸加工表面、可列印較大 面積的圖案(如使用陣列噴嘴系統)、可噴射多種流體材料、節省材料和少量 多樣製作等方面，故受到工業界青睞。由於目前工業的製程技術僅到微米 尺度，但隨著光電及其他產品的內部元件之微縮化的趨勢，因此微噴流列 印技術勢必需將尺度推展到奈米等級的噴流列印技術。微奈米噴嘴噴頭的 製造，可用雷射雕刻的方式，製作內部流道與噴嘴孔洞。 噴墨列印技術(或噴流加工技術)通常可分成兩類；連續式噴墨列印與供 需式 (Drop-on-demand) 噴墨列印(又稱為間歇式與脈衝式噴墨列印)。在連 續式噴墨系統中，壓力施加在噴嘴墨源上，使噴墨在離開噴嘴一定距離後 斷成墨滴，墨滴形成後再通過控制器使之感應帶電，由偏轉電極使帶電墨 滴偏轉，噴射到承印物(如基材)的指定位置，依照偏向方法的不同，可以分 為位元偏向式（Binary deflection）和多重偏向式（Multiple deflection）。其 缺點為解析度較低與僅適用於低黏度的材料，多用於印製大型海報與廣告 看材。

目前，噴流加工技術的主流是供需式噴墨列印系統，根據計算機指令 間歇地將墨滴噴到承印物上，形成圖案。依照驅動機構的不同，分為熱泡 式（Thermal Bubble）和壓電式（Piezoelectric）驅動兩大類。

圖 2-5 噴墨印表機技術之分類圖7 1. 熱泡式驅動技術 藉由一熱電阻(在幾微秒內即可迅速加熱達攝氏 350 到 450 度)，使水 氣化產生微氣泡，因體積膨脹而將原來空間中的墨水擠出於噴孔外。 此時加熱片又迅速冷卻使氣泡破裂，原擠出於噴孔外的墨水受到氣泡 破裂力量的牽引而形成墨滴，當墨滴噴出飛行至基材上時，便在基材上形 成墨點。而墨水則透過連通噴墨區與儲墨區的流道持續流入補充，每噴出 一個墨滴都是上述流程協同運作的結果。這種機構的噴墨速率每秒可達三

和側噴型（Side-Shooter）。 圖 2-6 熱泡式頂噴型噴墨技術示意圖7 2. 壓電式技術 為 EPSON 獨 立 開 發 的 驅 動 技 術 ， 其 關 鍵 在 於 晶 體 (Piezoelectric Material)的壓電特性(晶體通電時會產生膨脹的現象)，以多層壓電波穩定的 控制晶體上的電壓，使晶體準確產生膨脹將墨水噴出，所以墨點本身的產 生速率、大小都可以精準地控制。最後，當通電停止時便將墨水拉回7_。依 不 同 壓 電 產 生 的 形 變 機 制 ， 壓 電 式 又 可 分 為 收 縮 管 型 （ Squeeze Tube

Mode）、彎曲型（Bend Mode）、推擠型（Push Mode）及剪切型（Shear Mode）。

圖 2-7 壓電式推擠型噴墨示意圖7

由上面敘述可知，連續式與供需式噴流技術最大的不同在於液滴的帶 電與否；前者利用液滴帶電的特性，控制液滴飛行位置。後者則以精確的

定位噴嘴系統，控制液滴飛行位置。目前，商業化技術只到微米尺度。若 將來商業化技術擴展至奈米尺度。在奈米噴嘴噴頭的製造上8, 9_{，可用聚焦} 式離子束（Focused Ion-Beam）雕刻的方式，對矽材料製作內部的流道與噴 嘴的孔洞，噴嘴可以小於到 100 nm，藉由對離子束的調整與控制，可得到 理想的噴嘴幾何形狀。 隨著科技的進步，噴印製程也將朝向微縮化的方向發展，因此本研究 計畫將建立一奈米尺度之噴流模型進行研究。另外，在實驗研究中也確實 成功的在奈米尺度下進行噴印。

2.4. 研究動機與目的

目前與對奈米噴流加工技術相關的文獻僅止於學術上的研究，相關的 電腦模擬技術也離商業化軟體一大段距離10-17_{，其開發與應用的潛力受到各} 方的重視。因此，研究目標藉由分子動力學模擬開發奈米噴流加工製程的 軟體，建立下一世代之關鍵光電元件與圖案之設計基礎，並且建立奈米噴 流的流體動力學之理論分析能力。 噴墨列印技術為非接觸式的隔空置樣方式，由於未與基材表面接觸， 不會有針頭磨損、並可達高速均勻置樣之目的。另外，噴墨列印技術一般 均是單一噴頭負責單一溶液，所以並無清洗、乾燥流程，無交叉汙染之虞； 但也由於它不與基材接觸，所以在液滴移動期間會受到的環境影響變數較 多，需作較多的參數控制。 預定模擬的系統擬對供需式噴流技術進行研究，須以下列各點作為評 估的依據溶液從噴嘴噴出形成液滴的過程：需避免分子沾黏於噴嘴口表面 10

二次液滴斷裂 (Break-up)、各別液滴再團聚及液滴水平漂移。基本上，產生 上述五項現象都是研究奈米噴墨列印製程中要避免的情形。最佳化的設計 目標，是到達基材的液滴的大小一致，這些結果又會根據溶液分子本身的 交互作用力(與溶液的配方有關)、金屬原子與溶液分子的交互作用力(與噴 嘴材料的選擇有關)、加工條件與噴嘴幾何形狀等有關。 模擬上將假設系統處於真空狀態(目前，半導體製程也盡量使系統處於 真空狀態，可避免材料的污染)。以上的工作的完成，必須在知道各分子在 不同時間下的位置與速度等等數據，才能夠進行相關的分析，得到噴流的 最佳化的操作條件。上述的數據都可由分子動力學模擬的結果獲得，因此 分子動力學模擬與本研究的目的息息相關。

3. 文獻回顧

3.1. 水分子勢能模型

水是自然界最重要的溶劑，在化學及生物相關的領域也經常被使用18_， 因此水分子的性質及水分子對溶質之影響備受重視。 在分子動力學中，模擬結果之準確性主要根據是否能正確地描述分子 微觀結構及正確的統計性質的勢能模擬，對於極性甚強的水分子而言，其 氫鍵對於自身及溶質影響甚劇，因此提出數種不同的勢能模型。根據分子 點(sites)數、是否為柔性模型及是否包含極化作用(polarization effect)可分為 數類，常見的水分子模型有：

1. Three-sites model: SPC19, SPC/E20, SPC/Flexible18, 21, TIP3P22及 F3C23, 24

2. Multicenter model: TIP4P25及 TIP5P26

3. Polarizable model: SPC/FQ27, TIP4P/FQ27及 PLO128

每一種勢能的計算方式略有不同，如 SPC/E 與 TIP4P 一般為固定鍵長 與鍵角的演算法計算運動方程式，而 F3C 與 SPC/Fw 則使用具可撓性之鍵 長與鍵角18。然而實際上不論是哪一個勢能，都沒有辦法將水分子全部的性 質正確的計算出來，因此勢能的選擇將根據研究之系統的需求來進行選 擇。原則上考量勢能計算速度與性質計算之正確性做選擇。 對於大型的系統，如奈米噴流10-17_{、分子自組裝與生物方面} 29-35_之研究， 往往需要非常龐大的分子才能建立，所選擇之勢能對計算速度的影響甚 劇。因此分子動力學模擬在生物領域等較大之系統下 18_{，Three-sites model}

在 Three-site model 中，simple point-charge (SPC) model19由 Berendsen 等人所建立，是第一個能正確且以雙原子間勢能(pair potential)描述水溶液 之勢能。該水分子模型以三原子各別的質心 (atomic nuclei)計算分子間交互 作用力，而分子鍵長及鍵角透過固定(constrain)之不可撓曲模型(rigid model) 描述，凡得瓦爾力部分以常見的 Lennard-Jones 計算，並以 Columbic potential 計算分子間靜電力。

SPC 是相當優秀的水分子勢能模型，除了之前所提到的高計算速度之 外，多數的塊狀(bulk)水溶液性質也能獲得良好再現，求得與實驗值相近的

結果 36, 37 。然而十分重要的擴散係數(self-diffusion constant)及介電常數

(dielectric constant)並沒有辦法正確描述。

隨後 Berendsen 等人提出經過修正的 SPC model，extended SPC (SPC/E)

model20，該勢能與之前提出之 SPC model 僅有些許差異，相較於其它勢能

更能正確的描述所有水溶液的性質，原因歸咎於該勢能透過提高氫與氧之 部 分 電 荷 值 (partial charge) 修 正 水 分 子 自 身 的 極 化 能 量 (self-polarization energy)。然而這樣的修正雖然可以更加正確的描述水溶液，但是極化效應 (polarization effects)對於在水溶液周圍環境的物質影響甚劇，因此調整了部 分電荷值(partial charge)反而提高了在異相系統的不確定性，如薄膜與水油 兩相之系統。 相較於固定(constrain)分子模型，採用可撓曲之鍵長與鍵角的勢能演算 是較為合理的方式。故學者引進彈簧勢能籍以描述分子內作用力的行為， 並與過去提出之 SPC 與 TIP3P model 等結合。然而可撓曲(flexible)分子模型 所獲得的結果並沒有比非可撓曲(rigid)分子模型來得正確，反而增長計算時 間。因為在彈簧勢能的描述下，分子內相對運動行為變化的非常快，若想 描述該行為必需使用較小的時間步進(time-step)，因此，仍然有許多生物領 域的分子模擬依然採用固定(constrain)水分子模型。

其中還有一個非常優秀的勢能模型 F3C model，對於生物與工程之系統 可以說是絕佳的選擇。該勢能同樣採用可撓曲(flexible)鍵長、鍵角的演算 法，較合理的描述分子內行為。在非鍵結力的部分透過截斷法，並引進補 償 的 比 例 係 數 (A-scale) ， 減 少 凡 得 瓦 爾 與 庫 倫 靜 電 排 斥 能 量 (Repulsive energy)，補償因截斷半徑造成的吸引力損失。透過截斷勢能的方式描述長 程(long-range)勢能可大幅加快計算速度，節省計算時間。另一方面，在計 算結構性質(Structural properties)，動態性質(Dynamics properties)與勢力學性 質(Thermodynamics properties)都能準確地模擬計算而得，相較於其他勢能 僅得以針對某些特定性質的模擬計算，其描述水分子的能力優秀許多。

Levitt 等人23所發展的 ENCAD Package 之所以能夠準確的模擬各種性

質，是因為該研究團隊將龐大的分子模擬與實驗數據資料，藉由完整的理 論(量子計算與分子動力學)來建立了 F3C model。然而在近來的學術研究 中，該勢能模型仍多數運用在純水溶液特性之討論，如奈米液滴 38 等，幾 乎不見異相(heterogeneous)或溶液(solution)的研究，因此使用該勢能進行高 分子溶液方面的研究，將面臨缺乏可參考的可靠數據，對於以應用為目標 的研究相當的不利。 倘若能獲得一個水分子模型，可正確的計算大部分的性質，使用較合 理的可撓曲勢能並且是較為簡單的勢能保有較高的計算速度，將會是需要 大量水分子系統的最佳選擇，Yujie 等人18便提出這樣的勢能模型，該研究 團隊發現適當的調整水分子之鍵長與鍵角平衡參數，可使不正確的擴散係 數(self-diffusion constant)與介電常數性(static dielectric constant)質獲得改 善，同時不影響到原先其它性質的正確性。

衡長度有密切的關係，而介電常數與分子鍵角之平衡角度有關。當引進彈 簧勢能(harmonic potentials)，並提高平衡鍵長後，使原先 SPC 的徑向分佈函 數(radial-distribution function)向右偏移，並且產生較強的氫鍵(氫鍵有較長的 lifetimes)使水分子個別的移動性(mobility)下降，改善多數模型高估水分子動 力特性的結果。另一方面，提高平衡鍵角則會使介電常數數值下降。由於 調整平衡參數值並不會讓原先模型結構有太大的改變，因此除了之前提到 的兩個性質之外，依然能夠保有 SPC 在其他性質的正確性，由表 3.1 中可 以觀察到， SPC/Fw 除了某些性質的描述略遜於 SPC/E 之外，可以準確的 計算出所有的水分子特性，並且使用較合理的可撓曲勢能。 最為廣用的水勢能模型採用 Three-sites model，透過氫與氧各原子之質 心位置來計算交互作用力。在分子幾合結構上分為兩個 O-H 鍵長與一個 H-O-H 鍵角並給予適當的彈簧勢能函式，如下所示：

(

)

(

)

∑

⎪⎭ ⎪ ⎬ ⎫ ⎪⎩ ⎪ ⎨ ⎧ + ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ − ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ = − + − = _{∠ ∠} allpairs j i ij j i ij ij ij ij ij er HOH HOH bend OH OH bond ra r q q r r U k r r k U , 6 12 int 2 0 2 0 int 4 2 2 σ σ ε θ ( 3.1 ) 其中 ra

Uint 描述分子內勢能，Uinter則為分子間交互作用，r_OH0 與θ_∠0_HOH分別

為平衡鍵長與平衡鍵角，rij為 i 粒子與 j 粒子間的相對距離，εij與σij為

表 3.1 數種水分子勢能計算塊狀水溶液之性質。粗體字為最接近實驗數

3.2. 高分子與高分子溶液

由於噴墨列印技術的應用範圍相當廣泛，故選用應用範圍較廣的高分 子材料，除此分子動力學模擬之奈米噴流研究需要相當大的計算量，即便 電腦相關技術持續發展運算速度也不斷的提升，仍然不能滿足我們的需 求。因此本研究選用水溶性高分子聚氧乙烯作為研究材料。該材料可應用 金屬加工電鍍、電焊及油墨基質等，運用範圍甚廣，而且結構簡單同時具 有多種衍生材料，作為奈米高分子溶液噴流初步研究之材料相當的合適。

高分子材料聚氧乙烯 (聚環氧乙烷, Polyethylene oxide, PEO or poly oxyethylene, POE) 或聚乙二醇 (Polyethylene glycol, PEG)，合成上 PEO 以 環氧乙烷將碳、氧鍵結打斷開環，然後逐步形成高分子。PEO 與 PEG 在化 學結構上是完全相同的，差別在於 PEG 分子量在 20000 g/mol 以下，而 PEO 的分子量較大，因此兩者在物理性質方面有些許的差異，例如熔點與黏度， 而化學性質基本上是相同的。 聚氧乙烯透過分子動力學模擬已經被廣泛的研究，包括該材料之物理 性質及化學性質。除此之外，聚氧乙烯水溶液於表面吸附的特性也備受重 視，可應用在控制溶液中分子之聚集39，藉以保護蛋白質及細菌吸附在有機 或無機表面改善生物相容性(biocompatibility)40, 41_{。由於生物相關之研究需求} 大量的粒子及運算量，因此聚氧乙烯出現了各種不同尺度(length)模型的研 究。 分子動力學中，能夠正確的描述高分子材料的模型，熟知的主要有三 種 ， 分 別 為 Explicit atom (EA) model, United atom (UA) model 及 Coarse-grained (CG) model。EA model 為三者中描述最為詳細的模型，該模 型將各原子質心視為一個質點(site)並計算交互作用，包括分子鏈上所有的 氫原子都詳細描述。UA model 則為 EA model 的簡化模型，主要是將碳原

子及其鍵結之氫原子視為一個質點(site)，獲得較快速的計算效率。CG model 則是更加粗略描述分子鏈的模型，將數個 united atoms 或單體視為一分子群 並重新定義為一個質點(site)進行計算。在分子模擬中計算速度是相當重要 的考量之一，然而模型的選擇必需依研究之需求選定，倘若研究中希望得 知分子鍵中甲基的旋轉狀態等，那勢必要選擇 EA model，若希望能獲得較 為準確的性質計算速度較快則選擇 UA model。

對於聚氧乙烯而言，UA model 及 EA model 的研究較為透徹，而 CG

model42 相較於另外二者則顯的不足，即使擁有最適合本研究的高計算效

率。EA model 主要由 Smith43-46及 Neyertz47兩個研究團隊所建立，Smith 等

人所建立的勢能成功的求得聚氧乙烯的靜態及動態性質並與實驗獲得一致

的結果。Neyertz 等人則是重現結晶及熔融態的結構性質。隨後 Zon 等人48

以 Neyertz 所建立的勢能為基礎獲得 UA model 的勢能，成功計算散射函數 (scattering function)並與中子散射實驗結果一致。然而 Neyertz 及 Zon 等人 的研究中並沒有探討聚氧乙烯水溶液，也因此沒有 UA model 與水溶液間之 勢能的研究，因此，本研究選擇 Smith 等人所建立的聚氧乙烯勢能模型。

Smith 等人43於 1993 年以量子力學(ab inito)為基礎建立 1, 2-二甲氧基

乙烷的分子動力學參數，並利用結構或然率進行參數微調。1996 年 Smith 等人 44_{首次進行聚氧乙烯的模擬研究，並與小角度中子散射(SANS)實驗在} 各溫度下的結果進行比對，且發現到在結構與 1, 2-二甲氧基乙烷極為相 似。1998 年45_{以量子化學為基礎，利用 Lennard-Jones 描述聚氧乙烯與水之} 間的交互作用力。隨即於 1999 年46建立各種不同水勢能模型與聚氧乙烯之 間的作用力參數，並透過動態性質及分子結構比較彼此間的差異。

3.3. 奈米噴流

利用分子模擬進行液體噴流的研究方面，1993 年 Koplik 與 Banavar49 進行界面斷離（Rupture）的研究，圖 3-1 為他們研究液柱斷離形成液滴的 過程圖，系統由一開始的液柱逐漸斷離，然後慢慢調整，最後成為圓球； 他們改變不同的液柱半徑，並將液柱斷離的時間與線性理論做比較，得到 了不錯的結果。 圖 3-1 液柱斷離過程圖49 1998 年 Kawano50也進行了液柱在兩相斷離的研究，他加大了系統邊長 與液柱半徑的比值，所生成的液滴顆數並不像 Koplik 與 Banavar49 僅有單 顆，而下圖所生成的顆數較多，這是因為上圖的液柱半徑較大；而 Kawano 也將結果與 Rayleigh 的線性理論相比較，且得到相當不錯的結果。

圖 3-2 液柱斷離過程圖（Kawano）50

2000 年 Moseler 和 Landman10進行奈米噴流的研究，他們成功的模擬連

續噴流現象，進而討論整個液流的形成、穩定性及分離的現象，並且利用 連續力學的概念，將 LE（Lubrication Equation）加上一修正項項成為 SLE

（Stochastic Lubrication Equation），且由所得的 SLE 與分子模擬比較，得到

與分子模擬相當吻合的結果。圖 3-3 Landman 等人所模擬的奈米噴流情 形，這是考慮分子與噴嘴有潤溼（Wetting）現象，由圖中可看出液柱及液 滴皆出現非對稱情形，且伴隨著蒸發的分子；而在傳統連續力學模擬中， 一般常見的方式是將 Navier-Stokes 方程式簡化成一維軸對稱，因此在噴流 過程中有些現象無法忠實的呈現，雖然有人採用三維的模擬方式，可有效 的描述非對稱情況，但是當尺度小至奈米級時，對於分子間作用力所產生 的現象，仍然無法用巨觀的 Navier-Stokes 方程式來描述，圖 3-4 即為不同 方法所預測的情形，在液滴形成前，在液流頸部（Neck）MD 是呈現雙錐 形，利用連續力學的 LE 所得到的結果呈現細長的液流，SLE 模擬所得到的

圖 3-3 奈米噴流過程圖10

圖 3-4 液體噴流預測情形 10

2001 年 Goto51等人經由實驗與模擬進行研究，他們主要是利用雷射熔

射法（Laser Ablation），將一種有機分子注入另一種有機分子表面，如何成 功的將分子射出。

圖 3-5 雷射熔射法51 而在奈米噴流的實驗上，Voigt8, 9, 11, 52, 53等人做了一系列的研究，他們 由製作的奈米噴嘴，進行不同條件的奈米級蝕刻研究，他們發現，當基板 與噴嘴的距離小於噴嘴的半徑時，蝕刻的寬度會趨近於噴嘴的直徑，亦即 可產生約 100 nm 的高解析度；而因為單一噴嘴蝕刻速度太慢，所以為了克 服這個問題，他們同時也發展了奈米噴嘴陣列，如圖 3-6 所示。

而奈米噴嘴製作的過程，利用聚焦式離子束（Focused Ion-Beam）在金 字塔型的矽後方鑽孔，而噴嘴的幾何主要受到離子束流動分佈（Profile）的 影響，因此經由適當的調整控制，可得到理想的噴嘴幾何。 方得華(Fung)等人於 2003 年13_{研究 Argon 的噴流塗佈情形，研究中觀} 察不同噴嘴孔徑及移動基材對噴流塗佈的影響，並發現噴嘴內部能量隨時 間呈現波狀變化的情形。2004 年 14_{研究溫度及噴嘴孔徑對噴流的影響，發} 現除了能量之外壓力也呈現波狀變化的情形，並指出當系統溫度較低時有 利於流體集中不至於產生噴濺的現像。同作者於 2006 年16_{考量傾斜角度之} 噴嘴於基材的塗佈情形，此外也針對不同的噴嘴移動速競進行討論。 圖 3-7 Argon 於基材之噴流塗佈情形及傾斜噴嘴之塗佈情形 2005 年 Shin 等人15透過二維分佈圖 (2D contour) 發現噴嘴內部接近出 口處存在一分界點，當流體通過該界線後因為相變化 (Phase transition) 導 致包括軸向速度分佈、局部密度結構及能量等皆有明顯的變化。此外透過 不同收縮角度改變噴嘴模具的設計，然而並不會明顯的影響噴流的液柱及 液滴。此外提供一項重要的訊息，液柱的斷裂的情形主要由溫度來決定， 當溫度過低時需要較長的時間液滴才會發生斷裂；反之則會較快發生斷

裂，並且得到較短的斷裂長度。

4. 研究方法

4.1. 分子動力學理論

分子動力學以牛頓力學與漢米頓力學計算分子間作用力，配合適當的 數值方法預測粒子在空間中運動的軌跡、動量，主要有以下三項假設： 1. 所有粒子皆遵守牛頓運動定律。 2. 粒子之交互作用滿足疊加原理。粒子間因為固有的特性因此會受到一種 或多種以上的交互作用，而彼此間的關係透過量子力學或其他散射實驗 測量方法所求得之勢能方程式 (potential function) 來描述，進而求得對 粒子的效應，而該效應允許累加在一起。 3. 所有粒子均視為球型，並以質心描述其空間位置。 建立在這三項假設之下，配合由學者所提出的勢能方程式，以漢米頓 力學求得力場的函式 (force field)，並以適當的數值演算建立分子動力學模 擬流程。

4.2. 勢能模型

勢能函數是分子動力學的關鍵部分，選擇適當與否將大大的影響未來 的研究。必需要根據模擬的系統、運算的效率及統計的性質等做為考量。 本研究所選用之溶液為高分子溶液，並使用金屬作為噴嘴，因此接下來將 介紹溶劑、高分子及金屬的勢能函數。

4.2.1. 水分子勢能模型

本研究選用水作為溶劑，採用 Yujie 等人 18所提出 SPC/Fw model。其

分子勢能函式為被廣泛使用的 three-sites water model，各勢能函式如下：

1. 分子鏈拉伸作用力(Bond-Stretching interaction)

(

₀

)

2 2 OH OH bond bond r r k U = − ( 4.1 ) 將分子內鍵結透過虎克定律，以假想的彈簧勢能描述分子鍵運動行 為。其中kbond為彈性係數， 0 OH r 為平衡鍵長。 2. 分子鏈彎曲作用力(Angle-Bending interaction)

(

₂

)

2 2 HOH HOH bend bend k U = θ_∠ −θ_∠ ( 4.2 ) 如同拉伸作用力，不同之處在於描述鍵角的簡諧運動。其中kbend為彈性 係數，θ∠HOH為平衡鍵角。

3. 凡得瓦爾力(Van der Waal force)

為描述分子間短距(short-range)吸引力，為偶極-偶極(偶極力)、偶極-極 化(感應力)及分散力的總合。勢能函式如下： ∑ _⎥⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ − ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ =allpairs j i ij ij er r r U , 6 12 int 4ε σ σ ( 4.3 ) 其中

ε

與σ 為 Lennard-Jones 參數。

4. 庫侖靜電力(Coulomb force) 真空中兩靜止點電荷之間的交互作用力，電荷同號相斥、異號相吸。 ∑ =allpairs j i ij j i er r q q U , 0 int 4 1 πε ( 4.4 ) 其中qi則為部分電荷(partial charge)參數。上述之勢能參數列於表 4.1 中。 4.2.2. 聚氧乙烯勢能 本研究選用之高分子為水溶性聚氧乙烯，採用 Smith 等人 43所提出的 EA model，在鍵長及鍵角與 SPC/Fw 完全相同，由於長鏈分子具有超過三顆 以上的原子，必需考慮面與面的雙面角，在此以餘弦形勢描述： ( ) ∑ = − = 3 0 cos 5 . 0 k k Torsion a k U φ ( 4.5 ) 其中ak為雙面角能量參數，k則是餘弦波函數的摺疊數 (fold)，φ為雙 面角之角度。 分子間作用力同樣考慮庫侖靜電力與凡得瓦爾力，值得注意的是凡得 瓦爾力是以 Lennard-Jones exp-6 描述，方程式如下所示： ( ) 6 exp ij ij ij ij ij LJ r C r B A U = − − ( 4.6 ) 其中Aij、Bij、Cij為勢能參數。上述之勢能參數列於表 4.2中。

4.2.3. 金屬勢能

本 研 究 使 用 金 原 子 作 為 模 具 材 料 ， Morse potential model 、 EAM

(Embedded-atom method)或 TB (Tight-binding model)54, 55是最常見的金屬原

子勢能(屬於原子對原子間作用力)。其中 EAM 與 TB 為多體勢能(Many body potential)。多體勢能在計算兩原子作用力值時，除了需要考慮兩原子 之間的相對距離(如同 LJ model 一般)，同時考慮此二原子與其他的原子的 相對距離，故稱之為多體勢能，這是因為其勢能中加入了描述電子雲密度 的數學函式，為了能更精確地描述金屬原子，但相較來說，需要耗費更多 的電腦資源與更長的計算時間，對於需要大量金原子建立的系統而言，並 不是很好的選擇，因此本研究採用 Morse 勢能函式並選擇 Doyama 等人 56 所研究的金原子參數。該形式為較簡單的二體勢能，函式如下所示：

(

)

(

)

(

(

)

)

[

exp 2 r r0 2exp r r0

]

D U_Morse = − α − − −α − ( 4.7 ) 其中D為勢能能量參數。上述之勢能參數列於表 4.3中。 4.2.4. 其他分子間勢能 對於上述的三種材料之間的交互作用力，能否正確的使用適當的勢能 及參數，對於本研究有十分重大的影響。聚氧乙烯之所以選擇計算量較大 的 EA model 主要是考量到與水之間交互作用的描述，Smith 等人46_便針對 這部分進行研究，並提供了與各種水分子模型的參數，其勢能與參數列於表 4.4。同樣的，金原子與高分子溶液之間的交互作用力也十分重要，Dou 等

(

)

(

)

(

(

)

)

[

1 1

]

0 exp 2 O e 2exp O e O Au D r r r r U ₋ = − β − − −β − ( 4.8 )

(

)

(

2

)

0 exp 2 H e H Au D r r U ₋ =γ − β − ( 4.9 ) 該函式由 Spohr58所提出，與 Morse 勢能相似。考量到計算效率，因此 研究中會使用到截斷勢能，並使用開關函數(Switch function)59_{進行修正，讓} 勢能函數能平緩的衰退並於截斷距離歸零，該函式如下所示： ( ) ( ) ( )

(

) (

)

(

)

⎪ ⎪ ⎩ ⎪⎪ ⎨ ⎧ ≤ < − − + − ≤ = = ₋ − off on on off on off off on O Au O Au r r r r r r r r r r r r r S r U r S U 3 2 2 2 2 2 2 2 2 3 2 1 ( 4.10 ) 其中S( )r 為開關修正式，r_on為開啟半徑，r_off 為關閉半徑。 另一方面，金原子與聚氧乙烯的分子間作用力採用 Lorentz-Berthelot Mixing Rule 獲得，關係式如下所示：

(

A B

)

/2 AB σ σ σ = + ( 4.11 ) B A AB ε ε ε = ( 4.12 ) 使用 Lennard-Jones 12-6 勢能描述，參數列於表 4.3中。

表 4.1 SPC/Fw model 勢能參數表18 ( 0 )2 2 OH OH bond bond r r k U = − Bonds

(

2

)

/ /mol A kcal

k_bond & r_OH0

( )

A&

O-H 1059.162 1.012

(

₂

)

2 2 HOH HOH bend bend k U = θ_∠ −θ_∠ Bends

(

2

)

/ /mol rad kcal kbend

( )

deg 0 HOH θ H-O-H 75.90 113.24 ij j i ij ij er r q q r r U 0 6 12 int 4 1 4 πε σ σ ε + ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ − ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ = Nonbonded interaction

₍

₎

mol kcal/ ε σ

( )

A& qi

( )

e.c. O 0.1554253 3.165492 -0.82 H 0 -- 0.41 OH 0 -- --

表 4.2 Poly (ethylene oxide) 與 1,2-dimethoxyethane 勢能參數表42, 43 ( )2 0 2 r r k U b bond = − Bonds

(

2

)

/ /mol A kJ

k_bond & rij

( )

A&

0 C-C 2587.4 1.513 C-O 3094.0 1.390 C-H 2742.4 1.090 ( )2 0 2 θ θ θ ₋ =k U_bend Bends

(

kJ mol rad

)

kbend / / ijk

( )

A&

0 θ C-C-H 359.2 109.49 H-C-H 322.6 108.30 O-C-C 498.2 109.04 O-C-H 468.9 111.07 C-O-C 623.8 111.56 ( ) ∑ = − = 3 0 cos 5 . 0 k k Torsion a k U φ Torsions

(

kJ/mol

)

0 a a₁ a₂ a3 O-C-C-H 1.164 0.000 0.000 -1.164 H-C-O-H 1.164 0.000 0.000 -1.164 C-O-C-H 3.382 0.000 0.000 -3.382 O-C-C-O 10.886 0.209 -10.676 0.000 C-O-C-C 8.457 -4.187 -2.931 -1.340 ( ) ij j i ij ij ij ij ij LJ r q q r C r B A U 0 6 4 1 exp πε + − − = Nonbonded interaction

₍

₎

mol kJ Aij / Bij

( )

A&

(

)

6 /molA kJ C_ij & q_i

( )

e.c. C-C 62701.5 3.090 2682.9 -0.163 O-O 317547.0 4.063 1670.0 -0.256 H-H 11093.3 3.70 114.6 0.097 C-H 18087.0 3.415 578.8 -- C-O 141105.2 3.577 2116.7 -- O-H 59352.0 3.902 437.4 --

表 4.3 Au 相關勢能及參數56

(

)

(

)

(

(

)

)

[

exp 2 r r0 2exp r r0

]

D U_Morse = − α − − −α − Au

_{( )}

₋₁ A& α r &₀

( )

A D

( )

eV Au-Au 3.1132 2.8909 0.15479

(

)

(

)

(

(

)

)

[

1 1

]

0 exp 2 O e 2exp O e O Au D r r r r U ₋ = − β − − −β −

(

)

(

2

)

0exp 2 H e H Au D r r U ₋ =γ − β − H2O

( )

eV D0

( )

1 − A& β

( )

−1 A re & γ

( )

A& Au-H 6.68×10-2 1.3 3.7 0.2 Au-O 6.68×10-2 1.3 -1 -- ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ − ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ = 6 12 int 4 ij ij er r r U ε σ σ PEO

(

kJ /mol

)

ε σ

( )

A& Au-H 0.782 3.058 Au-C 2.436 3.280 Au-O 3.535 2.982 表 4.4 Polyethylene – Water 勢能參數表46 ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ − ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ = 6 12 int 4 ij ij er r r U ε σ σ PEO-SPC

(

kcal /mol

)

ε σ

( )

A& W-H 0.0750 3.0263 W-C 0.2335 3.2429 W-O 0.3389 2.9488

4.3. 初始化位置及速度

初始位置的排列並沒有一定的限制，不同的排列僅影響收斂至平衡的 速度，例如當對像為固體時，則依照物質在自然界的晶格排列方式給予初 始位置，若想要計算氣體的系統時，給予隨機的排列則會比較有利。一般 的作法，都是依照物質在固態時的晶格來進行排列，較常見的有四種堆積 結構方式，分別為簡單立方堆積 (SC)、體心立方堆積 (BCC)、面心立方堆 積 (HCC)及六方最密堆積 (HCC)。分子動力學模擬中，噴嘴多以單一金屬 建立，所以通常使用面心立方堆積排列。水的模型則是使用簡單立方堆積 或是面心立方堆積排列。另外，高分子的部分只要保持最穩定型態，分子 間取得適當距離即可。 建立初始位置後，藉由亂速給定每一顆粒子的速度，並將速度分量修 正到一個正負範圍之內，使粒子的速度呈現隨機的方向。透過以上的位置 及速度的配置，然後經過長時間的步進時間，即可使粒子混亂並呈現 Maxwell 速度分佈而擬似真實的狀態。 另一方面，在模擬上通常會將系統之質心速度歸零，使系統不會移動， 然而透過這種方式給予速度並不能確保系統總動量為零，因此必需在亂速 給定之後將總動量修正，其方程式如下： 1 / N new i i i k k i k m m m N ⎛ ⎞ =_⎜ − _⎟ ⎝

∑

⎠ v v v ( 4.13 ) 其中N為總粒子數，m為質量，v為粒子速度。

4.4. 系統類型與控制

在分子動力學中，會依照現實中的系統選定需要的系統類型，常見的 類 型 分 別 為 NVE ensemble 、 NVT ensemble 及 NPT ensemble 。 NVE ensemble 即粒子數 (N)、體積 (V) 及 內能 (E) 為恆定值，在奈米噴流的 系統中先進行 NVT ensemble 獲得噴流前系統，然後在不控溫的條件下進行 噴流模擬。 4.4.1. 系統體積控制 現階段的電腦計算效能不斷進步，但是仍然不足以提供分子模擬達到 巨 觀 的 粒 子 數 ， 在 有 限 的 計 算 量 之 下 經 常 採 用 週 期 性 邊 界 (Periodic boundary condition) 來模擬塊材 (bulk) 性質，如圖所示。透過該方法，僅 管無法求得正確值，但是這樣簡化的系統仍足以描述真實系統的特性。週 期性邊界的操作方式分為兩個部分： 1. 粒子位置的週期修正：當粒子超出系統時，必需將粒子修正會系統中維 持系統密度，方程式如下。 , 0, i i i i i i

L system boundary length when r L r r L when r r r L = > = − < = + ( 4.14 ) 2. 最小鏡像法：為了描述巨觀性質，當粒子間距離超過截斷半徑時，必需 將粒子間距修正到最短距離。目的是為了讓中心粒子周圍如充滿粒子的 巨觀系統，但實際上粒子數並沒有改變，所以當某一維度使用週期邊界 時，可以視該方向有無限延伸的複製系統所組成，若三維度皆使用，則

∑

∑

∑

∑

∑

∑

+ = < − = > 3 3 3 3 3 3 , 2 / , 2 / α α α α α α α α α α α α L r r  L r when L r r  L r when ij ij ij ij ij ij ( 4.15 ) , 2 x ij ij ij x L when x < − x =x +L 另外還有虛擬牆及實體牆邊界，在奈米噴流模擬中，噴嘴本身就是以 真實粒子組成的實體牆，而推動粒子噴出的驅動力通常也是來自真實推板 的推擠。 4.4.2. 系統控溫

1. Simple Velocity Scaling :

根據氣體動力論 d a d k a k B k T T E E T k E = = , , 2 3 ( 4.16 ) 其中Ta為實際溫度、Td為設定溫度。 a d i new i T T v v = ( 4.17 ) 該操作方法在科學上並不具有物理意義，依賴的是數學上的操作技 巧，因此該方法僅適用在不考慮模擬過程只求最終結果時使用。

2. Nose Hoover thermostat : 此方法修正牛頓運動方程式作進一步改寫，該方法如下 2 2 1 τ ζ ζ D B D B i i i i i i i T gk Q T gk m Q dt d  m = ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ − = − = =

∑

P p F p p v & ( 4.18 ) ζ 為摩擦係數，g為系統自由度，Q為溫控質量，τ 為特徵時間。在該 方法下，系統周圍如同存在一個恆溫裝置。當系統內的溫度過高時，則將 能量透過摩擦係數釋放能量，反之則提高能量，因此可以使系統在控溫過 程中維持能量恆定。

4.5. 數值演算法

在分子動力學中最常見的兩個演算法為 Gear’s 五階 predictor-corrector 及 velocity verlet 兩種數值方法，本研究使用 Gear’s 五階 predictor-corrector， 因為該方法可以獲得較精準的預測。演算的過程分成三個部分，分別為預 測 (Prediction)、計算 (Evaluation) 和修正 (Correction)，流程如下：

1. 預測 (Prediction) 預測的部分使用 Taylor-series 五階式進行演算，在初始化之後，藉由已 知的數值預測出下一個時間點的位置、速度及加速度等。方程式以矩陣方 式列於下方：

(

)

(

)

(

)

(

)

(

)

(

)

( )

( )

( )

_{( )}

( )

_{( )}

( )

_{( )}

( )

_{( )}

_⎥⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ ΔΔ Δ Δ Δ Δ Δ Δ Δ Δ Δ Δ Δ Δ Δ = ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ Δ + Δ + Δ + Δ + Δ +Δ + 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 5 4 3 2 2 1 ! 2 3 1 ! 3 2 1 ! 2 4 1 ! 4 3 1 ! 3 2 1 ! 2 5 1 ! 5 4 1 ! 4 3 1 ! 3 2 1 ! 2 1 0 0 0 0 0 1 0 0 0 0 1 0 0 0 1 0 0 1 0 1 r r r r v r r r r r v r Ⅴ Ⅳ Ⅲ Ⅱ ( 4.19 ) 2. 計算 (Evaluation) 藉由預測獲得的粒子位置配合勢能方程式計算出每顆粒子所受到的 力，並將作用力轉換成為加速度。 3. 修正 (Correction) 由於計算值的正確性要比預測所獲得的加速度可靠，因此將加速度修 正為較正確的值，並利用兩個加速度的差值配合適合的參數進行修正。

( ) ( ) ( ) ( )2 i, 2 e i, 2 i, 2 i, 2 -2! t r r r r R = Δ Δ Δ ≡ Δ ( 4.20 ) R r r r r v r r r r r v r Δ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ + ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ Δ Δ Δ Δ Δ = ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ Δ Δ Δ Δ Δ 5 4 3 2 1 0 5 1 ! 5 4 1 ! 4 3 1 ! 3 2 1 ! 2 5 1 ! 5 4 1 ! 4 3 1 ! 3 2 1 ! 2 a a a a a t t t t t t t t t t e α Ⅴ Ⅳ Ⅲ Ⅱ Ⅴ Ⅳ Ⅲ Ⅱ ( 4.21 ) 式中的a0、a1、a2、a3、a4及a5為 Gear’s 五階 Predictor-Corrector 的常 數值，會隨著演算法的階數 (q) 不同而改變，如下表所示：

表 4.5 三~五階 Gear Predict-Corrector Algorithm 參數

i α q=3 q=4 q=5 0 α 6 1 120 19 16 3 1 α 6 5 4 3 360 251 2 α 1 1 1 3 α 3 1 2 1 18 11 4 α -- 12 1 6 1 5 α -- -- 60 1

4.6. 加速計算方法

4.6.1. 牛頓第三運動定律 由於分子動力學假設分子運動符合牛頓運動定律，當使用作用力與反 作用力時，計算兩粒子間彼此的作用力時僅需要計算一次即可，兩粒子只 有受力方向不相同而已，因此可以減少一半的運算量。 4.6.2. 分子勢能截斷法 計算勢能以及粒子間作用力時，根據勢能方程式我們可以知道當粒子 間之相對距離大過一定的程度時，所獲得的數值非常小，因此，我們允許 設定截斷半徑 (Truncated potential)，當距離大過該半徑時將忽略該勢能數 值，換句話說，只要使用該方法，在超出範圍外的粒子可以不必計算，如 此一來我們將可節省大量的運算資源。

( )

r r     r r r r r U ij c ij ij ij ij > ≤ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ − ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ = , 0 , 4 6 12 σ σ ε _{( 4.22 )} 除此之外，使用截斷半徑除了節省粒子間交互作用的運算量，更重要 的是可以配合 Neighbor List 及 Cell-Linking 法，更加的提升運算效率。然 而，使用截斷半徑將引發其他的問題及限制條件，分為以下三點說明： 1. 作用力平移 (shifted-force potential)：當截斷勢能後，方程式將出現不連

續點。根據漢米頓力學，分子間作用力為勢能方程式微分並且取負值， 因此截斷勢能將會造成誤差，因此提出作用力平移修正不連續點，方程 式如下：

( )

⎪⎩ ⎪ ⎨ ⎧ > ≤ Δ + −∇ = c ij c ij ij r r   r r F U r , 0 , F ( 4.23 ) 其中ΔF =−F

( )

rc ，也就是 shifted-force。透過上式積分後可以獲得修正 後的勢能方程式。

∫

∫

∞ ⎟_⎠ ⎞ ⎜ ⎝ ⎛ + ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ − = s c c U u r r s dr dr dU dr dU dU 0 ( 4.24） ( ) ( ) ( ) [ ] ⎪ ⎩ ⎪ ⎨ ⎧ > ≤ ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ − − − = c c r c c s r r     r r dr dU r r r U r U r U c , 0 , ( 4.25）

圖 4-1 勢能 Lennard-Jones 12-6 shifted potential 示意圖

2. Long-range correction: 雖然在截斷半徑之外的勢能很小，但是若使用週 期邊界描述巨觀性質，則截斷半徑之外應該有非常龐大的粒子數被忽 略，僅管半徑之外勢能很小，但是若粒子數非常多時將會造成很大的誤

( ) ( )

∫

∞ = C r C _{U r g rr dr} N U 2 2πρ ( 4.26 ) 其中ρ為密度，g( )r 為徑向分佈函數於 4.8.3 中介紹。

3. Switch function: 當 使 用 截 斷 勢 能 時 ， 可 以 使 用 開 關 函 數 (Switch

function)59 進行較合理的修正。當分子間距離較靠近時，勢能不進行任 何的修正，至距離較遠處，勢能斜率變化較小時執行將勢能乘上一修正 式使勢能自然且平順的降低至零。以 4.2.4 提及的 Spohr 勢能為例，金 與氧之間的勢能衰退較慢，於 7 Angstrom 之前保持原勢能式，之後至 11Angstrom 之間開執行修正式，讓勢能函數能平緩的衰退並於截斷距 離歸零，如圖 4-2 所示：

圖 4-2 使用 Switch function 之 Spohr 勢能示意圖57

由圖 4-2 可以發現，Switch function 僅修正粒子相對距離較遠的勢能， 其勢能斜率變化小，經過修正後對整體影響並不大，因此該勢能優於作用 力平移法，其函式如下所示：

( ) ( ) ( )

(

) (

)

(

)

⎪ ⎪ ⎩ ⎪⎪ ⎨ ⎧ ≤ < − − + − ≤ = = ₋ − off on on off on off off on O Au O Au r r r r r r r r r r r r r S r U r S U 3 2 2 2 2 2 2 2 2 3 2 1 ( 4.27 ) 4.6.3. 鄰近列表法 分子動力學的計算量非常的龐大，僅管已經使用了截斷勢能法，運算 量隨著粒子數目的半方成正比關係，因此提出鄰近列表法 (Neighbor-list method) 加速計算。 鄰近列表法的基本觀念，在系統中的任一分子，其周圍的分子於短時 間內變動並不會太大，意即於此段時間內，此分子僅需和周圍鄰近固定的 分子計算作用力，可以不用對整個系統的分子作計算；當此方法使用在越 大的系統時，加速的效果更為顯著。 圖 4-3 鄰近列表法示意圖 鄰近列表法的使用方法，先將分子 i 列表半徑(List Distance ,r1)內的分

個時間間隔後再予與更新。 其中列表半徑和清單列表的頻率是相互影響的，選擇的決定依據在 於，位於列表半徑外的分子能否在下次清單更新前，進入影響分子的截斷 距離內，一般而言，當分子的移動速度越快，相對的列表更新頻率增高或 者是列表半徑要增加。 4.6.4. Cell-linking method 此方法使運算量僅與分子數目成一次方關係，大量減少運算量，並且 可與鄰近列表法做結合，更有效率的節省計算時所花費的時間。此方法的 概念是先將模擬的系統切割成若干等大小的方格(Cell)。 判斷出系統中每個分子所處的方格，當欲計算系統某一個分子受力情 形時，只需找尋在分子所屬方格周圍附近，未超過參考距離的數個方格內 的分子來判斷即可，不需要對整個系統的分子個別做搜尋判斷，而考慮使 用的距離應為截斷距離，但絕大部分此方法都會和鄰近列表法同時使用， 則參考距離就採用列表半徑。 圖 4-4 Cell-Linking 法示意圖

4.6.5. 列表勢能法 (tabulated potential method)

在系統進行模擬時，重復相似的迴圈運算之中，倘若遇到三角函數、 exponentiation 及 error function 時將耗費龐大的運算時間。本研究中水與高 分子溶液的結構特性研究當中必需使用 Ewald sum 計算庫侖力；高分子溶 液噴流研究中使用到 Morse 及 Lennard-Jones exp-6 勢能，大量的運用到這 類耗費時間的函數，若將這類函式或勢能和力直接列表，模擬時以查表方 式讀取數值，則可以省去很多的計算時間。然而列表的切割之粗細程度將

會影響數值的精確度，由文獻 60 中可知，當切割等分夠細時，實際上所產

生的誤差可忽略。本研究中包把 Erfc、Exp 及三角函數皆列表運算，並分割 為數萬等分， 圖 4-5 為 erfc 函數列表示意圖。

4.7. long-range interaction

勢能的運算中，包括 Lennard-Jones 及 Morse potential 皆屬於 short-range interaction，然而當勢能屬於 long-range interaction 時，一般的運算方式並不 能正確的描述，本節將討論如何處理 long-range interaction 的運算。 在 short-range interaction 計算塊狀性質，已經普遍的使用 4.4.1 提及的 週期性邊界條件，然而為了避免不正確的表面效應，因此允許使用 4.6.2 所 提到的截斷勢能法同時也帶來了提高運算效率的好處。然而對於 long-range interaction 而言，相同的處理方式將會發現到截斷空間之外的貢獻十分的 大，不合理的截斷將造成極大的誤差及不自然的分子運動。 當勢能允許使用截斷半徑時，都必需可以使用 long-range correlation 補 償截斷空間之外被忽略的貢獻。假設存在一勢能函式為 d r− 形式，並代入 4.6.2 的 long-range correlation 的勢能修正式中並積分，當d ≥3時，修正值將 為無窮大，意即無法使截斷勢能。對於這類分子間相對距力增加時，勢能 之交互作用降低很慢者視為 long-range interaction。因此，在 4.2 中提到的庫 侖力 (charge – charge, −1 r ) 的運算即屬於 long-range interaction。 若將庫侖力與凡得瓦力勢能相互比較，如圖 4-6 所示。在截斷距離處 可以看到 Lennard-Jones 幾乎為零，而庫侖力明顯存在相當大的數值，可想 而知，截斷之外的影響超乎我們的想像。當然，我們可以馬上想到一個直 接的處理方法，就是將系統直接放大，然而，可能放大至數百奈米之後可 能都無法符合要求。此外在不加入任何的加速方法時，計算時間與粒子數 之平方成正比，顯然這方法並不實際。