利用分子動力學研究不同操作條件下之高分子溶液的奈米噴流結構與行為

(1)

國立交通大學

應用化學系

碩士論文

利用分子動力學研究不同操作條件下之高

分子溶液的奈米噴流結構與行為

Study of the properties of nano-jet flow with

polymer solution under various operating

conditions via molecular dynamics simulations

研究生：蘇煜翔

指導教授：許千樹

指導教授：張榮語

(2)

摘要

軟性電子為一個新興的科技產品，近幾年來許多國家皆投入許多資金來進行研究，所潛在的商機可說是非常的龐大。於眾多技術中，噴墨技術於軟性電子的製程上佔有相當重要的地位。此外噴墨印刷技術也可以運用在電子元件組裝上，雖然目前仍無法跨入奈米尺度的印刷製程。倘若能達到奈米尺度的噴印，對未來的科技發展可說是一個相當大的突破。本篇研究希望藉由分子動力學模擬 (Molecular dynamics simulation) 來進行奈米尺度下的高分子溶液噴流研究，並透過不同的操作條件來觀察對噴流系統的影響。材料部份則使用結構較簡單、且應用範圍相當廣泛的水溶性高分子聚氧乙烯 (poly ethylene-oxide)。金屬勢能則選擇 Morse potential 來使用。水分子的部分使用 SPC/Fw model 來進行模擬。研究的過程分為三個階段： 1. 配合文獻資料，分別建立高分子溶液與金屬模型，驗證模擬之程式的正確性。 2. 將金屬建立成為噴嘴形狀，並以高分子溶液做為噴流材料，建立噴流系統。 3. 改變不同的操作變因來進行噴流模擬，並加以討論不同條件下對噴流的影響。

(3)

Abstract

Flexible electronics is a brand new technical product in these years. Due to the huge potential business opportunities, many countries have devoted many finances to develop this project. Among these technologies, micro scale jet printing technology plays an important role in the flexible electronics progress. In addition, the printing process can also be applied to assembling electronic components. Although the jet process can’t reach nanometer scale nowadays, once it can achieve one day, it will become a big breakthrough for the future technology development. In this research, we hope to concentrate on the nano-jet printing with a polymer solution as a printing material using molecular dynamics simulation by varying the operation conditions. Due to its simple structure and wide application, the water-soluble polymer, polyethylene oxide, was used by us as a printing material. Morse potential was chose to represent a potential model for the metal, and SPC/Fw model was chose to simulate the water.

The process of this study was divided into three steps：

1. Establishing the polymer solution and metal model, and then identifying the accuracy of our simulation program with the reference.

2. Using the metal to construct a nozzle and polymer solution as a printing material to establish the nano-jet system.

3. Discussing the different effects of the nano-jet system by changing the operation conditions.

(4)

誌謝

於碩士班的這兩年時間，首先是感謝許千樹老師與吳建興老師兩人的指導，使學生獲得研究所須的專業知識。同時也需感謝清華大學化工所的張榮語教授的共同指導，以及實驗室學長姐於程式上及觀念上的教導；其中最要感謝已畢業的洪健超學長當年花了相當多的時間在交接程式上，使這份碩士論文能順利的完成。雖然大部份的時間都在清大寫程式，於交大並沒有太大的交集。但也要感謝於交大的同學給與我的幫助，告訴我什麼時候要選課，什麼時候要注冊等。若沒有他們的幫助，或許我真的會忘了繳費而畢不了業。此外也要感謝當初於中正一起進來交大讀碩班的同學，雖然是不同的實驗室，但於無聊的時候可以去找他串門子，還可以順便去看超正的行政助理，讓自己有更好的心情來作研究；有時還可觀賞到他們實驗室最出名的肥皂劇，讓我於交大的生活不至於太過無趣。最後，也要感謝我的父母親對我的支持，讓我這段時間能心無旁騖的專心於研究上，也注接下來的學弟們於研究上能一切順利。

(5)

目錄

摘要... I Abstract ...II 誌謝... III 目錄... IV 圖目錄... VI 表目錄...X 符號說明... XI 1. 緒論 ...1 1.1. 研究目的與動機 ...1 1.2. 奈米科技產業發展概況 ...2 1.3. 軟性電子簡介 ...5 1.4. 噴墨印刷製程簡介 ... 10 2. 文獻回顧 ... 13 2.1. 奈米噴流 ... 13 2.2. 水分子勢能模型 ... 17 2.3. 高分子與高分子溶液 ... 19 3. 研究方法 ... 22 3.1. 分子動力學理論 ... 22 3.2. 勢能函數 ... 22 3.2.1.水分子勢能模型 ... 23 3.2.2.聚氧乙烯勢能 ... 24 3.2.3.金屬勢能 ... 25 3.2.4.其他分子間勢能 ... 25 3.3. 模擬流程 ... 29 3.4. 初始化位置及速度 ... 30 3.5. 系統類型與控制 ... 31

(6)

3.5.1.系統體積控制 ... 31 3.5.2.系統控溫 ... 32 3.6. 數值演算法 ... 33 3.7. 加速計算方法 ... 35 3.7.1.牛頓第三運動定律 ... 35 3.7.2.分子勢能截斷法 ... 35 3.7.3.鄰近列表法 ... 37 3.7.4.Cell-linking method ... 38

3.7.5.列表勢能法 (tabulated potential method)... 39

3.8. long-range interaction... 40 3.8.1.Ewald sum... 41 3.9. 性質計算 ... 42 3.9.1.溫度 ... 42 3.9.2.壓力 ... 43 3.9.3.徑向分佈函數 ... 43 4. 模擬系統 ... 45 4.1. 水與聚氧乙烯水溶液塊材系統 ... 45 4.2. 高分子溶液噴流系統 ... 47 5. 結果與討論 ... 49 5.1.SPC/Fw勢能模型之分子結構 ... 49 5.2.高分子溶液結構性質 ... 52 5.3.高分子溶液噴流 ... 55 5.3.1.不同推板速度對噴流系統的影響 ... 56 5.3.2.不同推板推擠方式對噴流系統的影響 ... 72 5.3.3.不同聚氧乙烯聚合度對噴流系統的影響 ... 87 6. 結論與未來展望 ... 96 參考文獻... 98

(7)

圖目錄

圖 1-1 軟性電子紙的優勢分析8_...7 圖 1-2 不同彎曲程度的軟性電子定義9_...8 圖 1-3 軟性顯示器產品發展進程9_...8 圖 1-4 EPSON 研發的軟性電子紙10_...9 圖 1-5 富士通開發之軟性電子紙及商品型錄展示11_...9 圖 1-6 噴墨印表機技術之分類圖13_...11 圖 1-7 微壓電式推擠型噴墨示意圖12_...12 圖 1-8 熱泡式頂噴型噴墨技術示意圖12_...12 圖 2-1 丙烷噴流模擬情形14_...13 圖 2-2 噴流之雙型結構14_...14 圖 2-3 不同噴嘴溫度造成液柱及密度中心偏移的現象49_...14 圖 2-4 不同形狀之噴嘴示意圖17_...15 圖 2-5 於 180K下之ARGON噴流情形17_...15 圖 2-6 ARGON之噴流塗佈於移動基材情形19_...16 圖 2-7 傾斜噴嘴之塗佈情形21_...16 圖 2-8 數種水分子勢能計算塊狀水溶液之性質表22_...19 圖 2-9 PEO應用示意圖...20 圖 2-10 聚氧乙烯於 300K下呈現縲旋結構38_...21 圖 2-11 PEO、WATER、LII不同混合方式所產生擴散係數40_...21 圖 3-1 分子動力學模擬流程圖...29

圖 3-2 勢能LENNARD-JONES 12-6 SHIFTED POTENTIAL示意圖...37

圖 3-3 鄰近列表法示意圖...38

圖 3-4 CELL-LINKING 法示意圖...39

圖 3-5 比較庫侖力與LENNARD-JONES勢能降低至零的程度...41

(8)

圖 4-2 DME水溶液平衡示意圖...46

圖 4-3 嘴模型幾何示意圖...48

圖 5-1 水分子 216 顆，溫度 298.16 K密度為1g/m3，SPC/FW MODEL徑向分佈函數之模擬結果圖。...50

圖 5-2 文獻22_{中比較實驗結果與SPC/FW MODEL的徑向分佈函數。}_...51

圖 5-3 DME水溶液於 318 K，DME莫耳分率(X_DME)為 0.18。分別為ODME-OWATER及ODME-HWATER 的徑向分佈函數。其中ODME-HWATER 的部分補償 2 的數值。...53

圖 5-4 文獻36_{中DME水溶液O} DME-OWATER及ODME-HWATER的徑向分佈函數。其中ODME-HWATER的部分補償 2 的數值。...53

圖 5-5 DME水溶液於 318 K。分別為DME莫耳分率(X_DME)0.18 及 0.72 的局部(LOCAL)與塊狀 (BULK)水溶液密度比值...54 圖 5-6 水獻36_{中，DME莫耳分率(X} DME)為 0.18 及 0.72 的局部(LOCAL)與塊狀(BULK)水溶液密度比值...54 圖 5-7 噴嘴與高分子溶液平衡後的切面示意圖...55 圖 5-8 噴嘴溫度 350 K，推擠速度 40 M/S之噴流瞬時圖...58 圖 5-9 噴嘴溫度 350 K，推擠速度 80 M/S之噴流瞬時圖...59 圖 5-10 噴嘴溫度 350 K，推擠速度 120 M/S之噴流瞬時圖...60 圖 5-11 推板推至相同位置(STEP 4)，不同推板速度之二維密度分佈圖...63 圖 5-12 推板推至相同位置(STEP 4)，不同推板速度之線性密度分佈圖...63 圖 5-13 推板推至相同位置(STEP 6)，不同推板速度的二維密度分佈圖...64 圖 5-14 推板推至相同位置(STEP 6)，不同推板速度之線性密度分佈圖...64 圖 5-15 推板推至相同位置( STEP 8)，不同推板速度之二維密度分佈圖...65 圖 5-16 推板推至相同位置(STEP 8)，不同推板速之線性密度分佈圖...65 圖 5-17 推板推至相同位置( STEP 4)，不同推板速度之二維壓力分佈圖...66 圖 5-18 推板推至相同位置(STEP 4)，不同推板速之線性壓力分佈圖...66 圖 5-19 推板推至相同位置( STEP 6)，不同推板速度之二維壓力分佈圖...67 圖 5-20 推板推至相同位置(STEP 6)，不同推板速之線性壓力分佈圖...67

(9)

圖 5-21 推板推至相同位置( STEP 8)，不同推板速度之二維壓力分佈圖...68 圖 5-22 推板推至相同位置(STEP 8)，不同推板速之線性壓力分佈圖...68 圖 5-23 推擠速度 40 M/S於溶液用盡後液滴自然飛行瞬時圖...70 圖 5-24 推擠速度 80 M/S於溶液用盡後液滴自然飛行瞬時圖...70 圖 5-25 推擠速度 120 M/S於溶液用盡後液滴自然飛行瞬時圖...71 圖 5-26 不同推擠速度下的液滴自然飛行瞬時圖...71 圖 5-27 CASE 1 條件下，推擠速度 40 M/S之噴流瞬時圖...74 圖 5-30 CASE 2 條件下，推擠速度 40 M/S之噴流瞬時圖...77 圖 5-31 CASE 2 條件下，推擠速度 80 M/S之噴流瞬時圖...77 圖 5-32 CASE 2 條件下，推擠速度 120 M/S之噴流瞬時圖...78 圖 5-33 CASE 1 條件下，於STEP 8 時不同推板速度之二維密度分佈圖...80 圖 5-34 CASE 1 條件下，於STEP 8 時不同推板速度之線性密度分佈圖...80 圖 5-35 CASE 2 條件下，於STEP 8 時不同推板速度之二維密度分佈圖...81 圖 5-36 CASE 2 條件下，於STEP 8 時不同推板速度之線性密度分佈圖...81 圖 5-37 CASE 1 條件下，於STEP 8 時不同推板速度之二維壓力分佈圖...82 圖 5-38 CASE 1 條件下，於STEP 8 時不同推板速度之線性壓力分佈圖...82 圖 5-39 CASE 2 條件下，於STEP 8 時不同推板速度之二維壓力分佈圖...83 圖 5-40 CASE 2 條件下，於STEP 8 時不同推板速度之線性壓力分佈圖...83 圖 5-41 CASE 1 條件下，不同推擠速度下的液滴自然飛行瞬時圖...85 圖 5-42 CASE 2 條件下，不同推擠速度下的液滴自然飛行瞬時圖...85 圖 5-43 不同推板條件下，推板速度與液滴自然飛行距離之線性分佈圖，其中CASE 0 為推板推至底的情形...86 圖 5-44 聚合度為 18 的聚氧乙烯水溶液，平衡後的切面示意圖...87 圖 5-45 聚合度為 18 的情形下，推擠速度 40 M/S之噴流瞬時圖...89 圖 5-46 聚合度為 18 的情形下，推擠速度 80 M/S之噴流瞬時圖...90 圖 5-47 聚合度為 18 的情形下，推擠速度 120 M/S之噴流瞬時圖...91 圖 5-48 聚合度為 18 時於STEP 8 時不同推板速度之二維密度分佈圖...92

(10)

圖 5-49 聚合度為 18 時於STEP 8 時不同推板速度之線性密度分佈圖...92 圖 5-50 聚合度為 18 時於STEP 8 時不同推板速度之二維壓力分佈圖...93 圖 5-51 聚合度為 18 時於STEP 8 時不同推板速度之線性壓力分佈圖...93 圖 5-52 聚合度為 18 時不同推擠速度下的液滴自然飛行瞬時圖...94 圖 5-53 CASE 1 條件下，聚合度為 18 時不同推擠速度下的液滴自然飛行瞬時圖...94 圖 5-54 CASE 2 條件下，聚合度為 18 時不同推擠速度下的液滴自然飛行瞬時圖...95 圖 5-55 所有不同條件下，推板速度與液滴自然飛行距離之線性分佈圖，其中N代表聚氧乙烯的聚合度...95

(11)

表目錄

表 1-1 全球奈米科技研發、中央政府投入經費3_...3

表 1-2 1976~2002 年奈米專利人權國家專利數目2_...4

表 1-3 2007～2017 年間的軟性顯示器（ FLEXIBLE DISPLAY ）市場預測7_...6

表 3.1 SPC/FW MODEL勢能參數表(YUJIE, HARALD ET AL. 2006)...26

表 3.2 POLY (ETHYLENE OXIDE) 與 1,2-DIMETHOXYETHANE勢能參數表(SMITH, JAFFE ET AL. 1993; CHUNXIA, PRAVEEN ET AL. 2006)...27

表 3.3 AU相關勢能及參數(IMAFUKU, SASAJIMA ET AL. 1986)...28

表 3.4 POLYETHYLENE – WATER勢能參數表36_...28

表 3.5 三~五階GEAR PREDICT-CORRECTOR ALGORITHM 參數...35

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

(12)

符號說明

bond U ：分子鏈拉伸勢能 bend U ：分子鏈彎曲勢能 OH r ：水分子鍵長 0 OH r _{：水分子平衡鍵長} HOH ∠ θ _{：水分子鍵角} 0 HOH ∠ θ _{：水分子平衡鍵角} er Uint _{：分子間勢能} ra Uint _{：分子內勢能} σ ε、 _{：Lennard-Jones 12-6 勢能參數} i q _{：部分電荷大小 (partial charge)} k a _{：雙面角勢能參數} φ_{：雙面角角度} ij ij ij B C A、、： Lennard-Jones exp-6 勢能參數 D：Morse potential 勢能參數 ( )r S _{：開關函數(Switch function) 勢能} 修正式 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 _：動能

(13)

1. 緒論

1.1. 研究目的與動機

物質在奈米尺寸下，會有別於原有尺度下的物理、化學或生物特性。奈米科技便是運用這方面的知識，在奈米尺寸的微小世界中操作、控制原子或分子組合成新的結構，以便展現新的機能與特性。並以其為基礎設計、製作、新的材料和器具，使它產生全新功能，並加以利用的技術總稱。奈米科技涵蓋的領域甚廣，從基礎科學橫跨至應用科學，包括物理、化學、材料、光電、生物及醫藥等1_。目前全球的工業發展皆趨向於奈米科技加工製程技術，奈米加工製程將取代現今的加工技術。然而，許多於巨觀尺度下透過連續力學所觀察到的特性於奈米尺度下確會出現截然不同結果。對於奈米尺度之下的知識與技術，仍然有許多不了解的問題及尚未發覺的現象。此時電腦模擬便成為一種相當方便的方法。分子動力學 (Molecular dynamics simulation)以微觀的角度來進行研究分析，於奈米尺度的研究上便扮演著一個相當重要的角色。目前已被運用在如高分子之流變學、奈米加工製程及生物工程等多種領域的研究。利用分子動力學進行電腦模擬，其具有以下優點： 1. 實際的實驗中往往需要價格十分昂貴的器材，且實驗往往需要花費大量的時間、材料以及人力。對於產品的開發將會是個障礙，甚至可能延誤商機。電腦模擬具有相當低的成本，可以確實的降低實驗的成本。 2. 於實驗上要今精確的統計奈米尺度下的數據有一定的困難度，透

(14)

過模擬的方式不但不受限於儀器、設備，甚至能突破實驗無法做到的操作條件，彌補實驗不足的地方。 3. 近年來電腦計算速度與記憶體容量大幅提昇，隨著電腦效能的提昇，分子動力學模擬的精度及速度也因此跟著大幅度的提升。目前已有許多研究資料顯示，分子動力學所模擬出來的結果與實驗所觀察得到的現象都能高度吻合4 , 5，即便不能獲得精確值，但在定性方面同樣能獲得良好的參考。噴墨技術於 1970 年便開始有相關的研究，至今已是一個相當成熟的技術。由於噴墨過程中不會與基材表面接觸，因此不會有針頭磨損的情形發生。且低成本、快速、高精確的定量微液滴體積、可精確的將液滴控制於所需的定位等優點，提供了一個現今工業相當具潛力的材料微分配技術。目前的噴墨列印系統是屬於微米等級，奈米噴流加工的僅止於學術上的研究，但其開發與應用的潛力均受到各方的重視。本研究預計使用分子動力學模擬來探討奈米噴流過程，並藉由改變不同操作變因來探討對噴流系統造成的影響。

1.2. 奈米科技產業發展概況

隨著科技的進步，奈米技術的知識已經建立了許多的應用技術，例如奈米結晶材料、奈米粉體、複合材料、塗佈及奈米纜線的相關技術，運用這些技術研發出許多應用產品。奈米科技對外來將帶來巨大的影響力，各國都全力投入該科技之研發工作。綜觀國際目前奈米科技發展趨勢，美國在奈米結構與自組裝技術、奈米粉體、奈米管、奈米電子元件及奈米生物技術有顯著發展；德國則在奈米材料、奈米量測及奈米薄膜技術略有領先；而日本則在奈米電子元件、無機奈米材

(15)

料領域較具優勢2。奈米科技的重要性是全球有目共賭，世界各主要經濟體無不視奈米科技為國家科技發展之重要領域，並且投入各方面資源，以及各種策略規劃來提升其奈米科技之發展。其中又以美國與日本的投入最為積極，佔全球奈米研發資金 50%以上，由此也凸顯出美日兩國在奈米科技研發上仍居主導地位。自 1997 年起美國對於奈米科技發展投入最多，其次是歐盟與日本。以全球投資趨勢而言，至 2006 年為止，世界各主要經濟體投入奈米科技的資金相較 1997 年時已達 10 倍之多。台灣自 2003 年正式開始推動奈米國家型科技計畫，並規劃在 2003 年到 2008 年間投入約新台幣 230 億元經費於奈米科技發展；藉由「奈米國家型科技計畫」有效整合產學各界力量及資源，建立我國發展學術卓越和相關應用產業所需要之奈米平台技術3。表 1-1 全球奈米科技研發、中央政府投入經費3

(16)

過去廿年來，奈米科技專利的成長率驚人。奈米科技專利數量最多的前十名國家分別為美國、日本、法國、英國、台灣、韓國、荷蘭、瑞士、義大利和澳洲。美國、法國、日本、英國、瑞士、荷蘭和義大利自一九七零年開始發表奈米技術專利，到一九九零年初期，韓國和台灣陸續跟進。最近五年成長最快的領域則是化學和製藥業，其次便是半導體元件。整體來說，全球開始愈來愈注意奈米科技知識。表 1-2 1976~2002 年奈米專利人權國家專利數目2

(17)

1.3. 軟性電子簡介

光電元件中，尤其以軟性電子 (Flexible Electronics) 相當具有潛力。工業技術研究院也看好軟性電子產業潛力，投入大量技術研發，並與 50 餘家廠商共組軟性電子產業推動聯盟、軟性天線印製及系統研發聯盟、連續式軟性液晶薄膜研發聯盟，並共同發展軟性記憶型顯示器材料開發驗證計畫、軟性電子關鍵設備技術整合計畫。目前已獲 217 項關鍵專利，居世界第一。經濟部表示，未來 3 年為軟性電子最佳發展時機。工研院電子所指出，2009 年全球軟性電子產值達 58 億美元，2015 年全球產值估計可達 237 億美元。軟性電子技術可讓電子產品的應用領域進一步擴展，使電子產品具備輕薄短小、符合人因工程、便利穿戴、柔軟安全、可自由捲曲的特性，發展出適合行動的數位化多功能電子生活產品。目前全球經濟情勢不佳，今年台灣面板產業面臨嚴峻挑戰，軟性電子相關技術可立即與產業現有製程銜接，生產軟性顯示器，經濟部技術處長杜紫軍表示，介於紙與非紙之間的軟性電子核心技術尚在萌芽階段，每年產值平均複合成長率達 65%，具廣泛應用商機6。據 Displaybank 於 2007 年 5 月 10 日公佈 2007 - 2017 年軟性顯示器（ Flexible Display ）市場預測， 2008 年以前市面上的軟性顯示器為 7 英寸以下的電子紙，主要用於電子書籍，2008 年以後電子紙用途將逐漸擴大，突破應用產品的尺寸技術延伸至 30 英寸以上(大型看板)。2010 年電子紙市場規模將達到 2 億 8000 萬美元，預計 2013 年有可能出現 50 英寸以上的電子紙。另外，有機 EL 面板與液晶面板領域的軟性顯示器也將出現；市場規模方面，到 2015 年時將達到 59 億美元，到 2017 年時將擴大至 122 億美元7。

(18)

表 1-3 2007～2017 年間的軟性顯示器（ Flexible Display ）市場預測7 軟性電子是一種技術的通稱，是一種建置在薄塑膠片、薄玻璃或薄金屬片的軟性或可撓性基材上的微電子元件與材料技術。軟性電子使用非晶矽、低溫多晶矽以及有機材料作為半導體材料，由於材料的可塑性，因此可以透過噴墨列印的印製式製程將導電材料佈設在塑膠、薄型玻璃或薄型金屬等可撓性的基材上，讓電子元件具有可撓性以及大面積製作產出。此外軟性顯示器在螢幕畫面不變時，是不耗電的，只有在更換螢幕畫面時，才需要耗電，與目前平面顯示器需要一直耗電的狀況是不相同的，也正因為具有此特性，使得軟性顯示器的耗電量非常低，進而也減少了電池的更換頻率，且觀看角度不受限制。由此可見，奈米噴流列印技術未來將在軟性電子製程技術中扮演極為重要的角色。

(19)

圖 1-1 軟性電子紙的優勢分析8 而在所有軟性電子應用領域當中，以軟性顯示器的研發進展與產品商業化也最為快速，投入廠商數目也較其他領域來的多；其主要原因是軟性顯示器的應用範圍相當廣泛，而由玻璃所構成的傳統平面顯示器在重量、厚度、固定型態等，在某些應用裝置的應用上造成一定的負擔，直接促使軟性顯示器的切入機會。除此之外，過去靜態資訊的顯示也可藉由軟性顯示器在可彎曲、成本低等特性，創造出許多新的產品應用，例如包裝外盒、價格標籤、藝術品創作等。

(20)

圖 1-2 不同彎曲程度的軟性電子定義9

(21)

圖 1-4 EPSON 研發的軟性電子紙10

(22)

1.4. 噴墨印刷製程簡介

噴印加工製程的優點在於能精準控制液滴的量與位置、生產速度快速、成本較其他製程低廉、列印圖案可具有高解度、減少半導體微影製程的次數、不接觸加工表面、可列印較大面積的圖案(如使用陣列噴嘴系統)、可噴射多種流體材料、節省材料和少量多樣製作等方面，故受到工業界青睞。由於目前工業的製程技術僅到微米尺度，但隨著光電及其他產品的內部元件之微縮化的趨勢，因此微噴流列印技術勢必需將尺度推展到奈米等級的噴流列印技術。噴墨列印技術通常可分成兩類；連續式(Continuous)與供需式 (Drop-on-demand) 。在連續式噴墨系統中，主要是利用壓力施加在噴嘴墨源上，使噴墨離開噴嘴。墨滴形成後再藉由靜電或電磁方式，使液滴感應帶電，噴射到指定位置。依照偏向方法的不同，可以分為位

元偏向式（Binary deflection）和多重偏向式（Multiple deflection）。其

優點為列印速度快且適用的化學藥品多；缺點為解析度較低且僅適用於低黏度的材料，目前多用於印製大型海報與廣告看材。目前，噴墨列印技術的主流是供需式噴墨列印系統，根據計算機指令間歇地將墨滴噴到承印物上，形成圖案。依照驅動機構的不同，可分為熱泡式（Thermal Bubble）和微壓電式（Piezoelectric）驅動兩大類。

(23)

圖 1-6 噴墨印表機技術之分類圖13 1. 壓電式技術為 EPSON 獨立開發的驅動技術，其關鍵在於晶體通電時會產生膨脹的現象。微壓電技術是利用壓電晶體在兩端電壓變化的情況下能以固定頻率震動的特性，透過電信號驅動噴墨室內的壓電晶體運動，壓迫墨水噴出。此技術的好處在於能夠精確控制墨滴的產生速率、大小、位置，列印精度很高；而且，使用這種技術的噴頭由於不需要加熱，因此受到腐蝕的機會較小，壽命要長一些。但是這種技術也有缺點存在，一是隨著列印頭壓電晶體數量的增多，其集成度不可避免的要越到一定的困難，而噴嘴數量的多少在一定程度下決定著列印速度的高低；而其另一個缺點就在於這種噴頭製造成本相對較高，所以導致了其印表機機身價格也較高一些。壓電式又可分為收縮管型

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

(24)

圖 1-7 微壓電式推擠型噴墨示意圖12 2. 熱泡式驅動技術熱泡式噴墨方式的工作原理與壓電式有著本質的區別，它是透過加熱電阻對噴頭加熱使墨水腔內的空氣受熱膨脹，產生壓力把墨水 “擠出＂噴嘴，形成噴向列印介質的墨滴。此種列印方式的好處在於實現比較簡單，因此成本較為低廉，有助於其機身售價的降低。而且現在隨著光學蝕刻工藝的發展，使用這種方式的印表機可以輕易的實現高密度的噴嘴排列，使得其噴嘴數量遠遠超過採用微壓電技術的產品，而更多的噴嘴就意味著更大的單程列印覆蓋率，也意味著更快的列印速度。而與其優勢相比，熱泡式噴墨方式的缺點同樣明顯。採用這種列印方式的產品在墨滴噴出時不能向微壓電技術那樣很好的控制墨滴的位置，因此在列印精度上要差一些，其次，由於噴嘴在列印時需要不斷地加熱，因此其更容易遭到腐蝕，這直接影響著列印頭的壽命。熱泡式主要可分為頂噴型（ Roof-Shooter ）和側噴型（Side-Shooter）12。圖 1-8 熱泡式頂噴型噴墨技術示意圖12

(25)

2. 文獻回顧

2.1. 奈米噴流

Landman14於 2000 年於 Science 發表奈米噴流之研究，模擬丙烷

在奈米尺度下的噴流情形。透過其模擬結果，可發現當液柱即將斷裂時產生了軸對稱的雙錐形形態，並將連續力學中 Lubrication equation(LE)修正為 stochastic lubrication equation (SLE)描述流體之自由表面。此外也分析了 wetting 對噴流的影響。

(26)

圖 2-2 噴流之雙型結構14

Landman 的學生 Wei Kang 於 2008 年49嘗試將噴嘴分為上下兩部

份，並分別設定為不同的噴嘴溫度。結果發現液柱與密度中心皆會偏移向噴嘴溫度較低的那一部份。圖 2-3 不同噴嘴溫度造成液柱及密度中心偏移的現象49 2005 年 Shin 等人 17研究不同的噴嘴結構及溫度對噴流的影響。不同的噴嘴結構對於噴流並不會有太大的影響。另外，液柱的斷裂的情形主要由溫度來決定，當溶液溫度過低時液滴會較慢斷裂。反之，則會較快發生，並且會具有較短的斷裂長度。

(27)

圖 2-4 不同形狀之噴嘴示意圖17 圖 2-5 於 180K 下之 ARGON 噴流情形17 Fang, T.-H. 等人則對於 Argon 的噴流情形作了一系列的研究。先於 2003 年 19_{研究 Argon 的噴流塗佈情形，研究中使用不同噴嘴孔徑} 及移動基材對噴流塗佈的影響，並發現噴嘴內部能量會隨著時間呈現波狀變化的情形。

(28)

圖 2-6 ARGON 之噴流塗佈於移動基材情形19 2004 年 20研究溫度及噴嘴孔徑對噴流的影響，發現除了能量之外壓力也呈現波狀變化的情形，並指出當系統溫度較低時有利於流體集中不至於產生噴濺的現象。於 2006 年 21 _{則進一步討論傾斜角度之噴嘴對於基材的塗佈情} 形，此外也針對不同的噴嘴移動速度進行討論。圖 2-7 傾斜噴嘴之塗佈情形21

(29)

2.2. 水分子勢能模型

水是自然界中最廣泛被應用的溶劑，在化學及生物相關的領域經常被使用22_{，因此水分子的性質及水分子與溶質之間的交互影響受當} 相當程度的重視。對於極性甚強的水分子而言，氫鍵對於自身及溶質性質影響甚劇，模擬結果之準確性主要根據是否能正確地描述分子微觀結構及正確的統計性質的勢能模擬，因此有許多學者提出不同的勢能參數來描述水分子。每一種勢能的計算方式略有不同，如 SPC/E 與 TIP4P 使用固定鍵長與鍵角的方法計算運動方程式，而 SPC/Fw 則使用具可撓性之鍵長與鍵角 22_{。依據分子點(sites)數、是否為柔性模型及是否包} 含極化作用(polarization effect)可大致將水分子模型分為三類，其常見的水分子模型如下：

1. Three-sites model: SPC23, SPC/E24, SPC/Fw22 , 25, TIP3P26及 F3C27 , 28

2. Multicenter model: TIP4P29及 TIP5P30

3. Polarizable model: SPC/FQ31, TIP4P/FQ31及 PLO132

對於一些較大型的系統，如奈米噴流、奈米壓印或生物基因工程方面的研究，往往需要使用到相當龐大的分子量，勢能的計算速度往往是最需被考慮的。因此分子動力學在模擬較大的系統時 22_， Three-sites model 往往是最常被使用的勢能，該勢能相較於其它類型具有最佳的計算效率。不同的水分子模型對於一些特定的性質會有良好的穩合度，但也會有某些性質會產生較大的偏差。倘若能獲得一個水分子模型，可正

(30)

確的計算大部分的性質、使用較合理的可撓曲勢能並且是較為簡單的勢能保有較高的計算速度，將會是需要大量水分子系統的最佳選擇，

Yujie 等人22便提出這樣的勢能模型，稱之為 SPC/Fw moedl。該勢能

模型以早期提出的 SPC model 為基礎，同時引進可撓曲鍵長與鍵角勢能的觀念，並適當的調整水分子之鍵長與鍵角平衡參數，使擴散係數 (self-diffusion constant)與介電常數(static dielectric constant)性質獲得大幅度的改善，同時於其他性質上也有一定的正確性。由圖表 3.7 中可以觀察到，SPC/Fw 除了某些性質的描述略遜於 SPC/E 之外，均可以準確的表現出所有的水分子特性。同時於其研究中發現，水的擴散係數與平衡鍵長之長度有一定的關係，平衡鍵長增長會使得氫鍵交連 (network)的機會增加，使得水分子難以移動而造成擴散係數下降。介電常數則是與平衡鍵角大小有關，較大的平衡鍵角會使得分子偶極距 (dipole moment)變小，使得介電常數隨之變小。 SPC/Fw 水勢能模型也是採用 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 σ σ ε θ ( 2.1 ) 其中 ra Uint 描述分子內勢能，Uinter則為分子間交互作用，rOH0 與 0 HOH ∠ θ 分別為平衡鍵長與平衡鍵角，rij為 i 粒子與 j 粒子間的相對距離，εij與 ij

(31)

圖 2-8 數種水分子勢能計算塊狀水溶液之性質表22

2.3. 高分子與高分子溶液

聚氧乙烯分子式為 H-(-OCH2CH2-)n-OH，外觀為白色粒狀粉末。 Tg 點為 65℃ ~ 67℃。完全溶於水，並會溶於某些有機溶劑。其溶液在低濃度下具有很高的黏性。具有柔軟性、高強度的熱可塑性樹脂。可以用壓延、擠壓、鑄塑等方式進行加工耐細菌侵蝕，不會腐敗，在大氣中的吸濕性通常不大。該材料可應用於金屬加工電鍍、高分子加工及噴印墨水穩定劑等範圍。於 2.1 的文獻回顧中，可發現已有許多學者投入奈米噴流的研究中。但所使用的大部份都為單一種材料，極少人使用混合溶液來進行研究。聚氧乙烯具有結構簡單與具有親水性等特性，因此本研究選來進行高分子溶液噴流的材料。

(32)

圖 2-9 PEO 應用示意圖 Grant D. Smith 等人對聚氧化乙烯做了一系列相關的研究。首先於 1993 年 33_{進行分子動力學模擬的研究，透過量子力學方法的計算} (ab initio)來獲得聚氧乙烯的作用力場參數，同時與實驗值進行比較。在構形方面，過去的研究中 DME 的構形傾向 tgg > tgt >ttt，因此模擬也必需獲得一致的結果。1996 年 Smith 等人 34首次進行聚氧乙烯的模擬研究，並與小角度中子散射(SANS)實驗在各溫度下的結果進行比對。1998 年35以量子化學為基礎，利用 Lennard-Jones 描述聚氧乙烯與水之間的交互作用力。隨即於 1999 年建立各種不同水勢能模型與聚氧乙烯之間的作用力參數，並透過動態性質及分子結構比較彼此間的差異。於 1996 年 Tasaki 開始研究聚氧乙烯的溶劑效應，先後加入了水以及苯38 , 39進行討論，研究中發現於 300K 時聚氧乙烯在溶劑中會呈

(33)

現縲旋狀的結構。

圖 2-10 聚氧乙烯於 300K 下呈現縲旋結構38

在聚氧乙烯中加入 Li 離子可藉由 Li 離子的移動性而增加導電度，可用來做固態高分子電解液(Solid polymer electrolytes)。所以於 2007 年 Tao, Z.; Cummings, P. T.嘗試於 PEO/Li 離子中再加入水分子進

行研究40_{，發現加入水以後 Li 離子的擴散能力提昇了約 50 倍以上。}

能使導電度大幅度的上升。

(34)

3. 研究方法

3.1. 分子動力學理論

分子動力學是以牛頓力學與漢米爾頓力學計算分子間作用力，並配合適當的數值方法預測粒子在空間中運動的軌跡、動量，主要有以下三項假設： 1. 所有粒子皆遵守牛頓運動定律。 2. 粒子之交互作用滿足疊加原理。 3. 所有粒子均視為球型，並以質心來描述其空間位置。建立在這三項假設之下，配合由學者所提出的勢能方程式，以漢米爾頓力學求得作用力函式(force field)，並以適當的數值演算建立分子動力學模擬流程。

3.2. 勢能函數

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

(35)

3.2.1. 水分子勢能模型

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

model22。其分子勢能函式為被廣泛使用的 three-sites water model，各

勢能函式如下： 1. 分子鏈拉伸作用力(Bond-Stretching interaction)

(

₀

)

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

(

₂

)

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

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

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

(36)

4. 庫侖靜電力(Coulomb force) 真空中兩靜止點電荷之間的交互作用力，電荷同號相斥、異號相吸。 ∑ =allpairs j i ij j i er r q q U , 0 int 4 1 πε ( 3.4 ) 其中qi則為部分電荷(partial charge)參數。上述之勢能參數列於表 4.1 中。 3.2.2. 聚氧乙烯勢能 本研究選用之高分子為水溶性聚氧乙烯，採用 Smith 等人(Smith, Jaffe et al. 1993)所提出的 EA model，除了前面提到的鍵長及鍵角外，由於長鏈分子具有超過三顆以上的原子，必需考慮面與面的雙面角 (Torsion angles interaction)，在此以餘弦方式來描述：

0 2 ( ) 2 bond Bond ij ij k U = r −r ( 3.5 ) 0 2 ( ) 2 bend

Bend ijk ijk

k U = θ −θ ( 3.6 ) ( ) ∑ = − = 3 0 cos 5 . 0 k k Torsion a k U φ ( 3.7 ) 其中ak為雙面角能量參數，k則是餘弦波函數的摺疊數 (fold)，φ 為雙面角之角度。分子間作用力同樣考慮庫侖靜電力與凡得瓦爾力，值得注意的是凡得瓦爾力是以 Lennard-Jones exp-6 描述，方程式如下所示： ( ) 6 exp ij ij ij ij ij LJ r C r B A U = − − ( 3.8 )

(37)

其中Aij、Bij、Cij為勢能參數。上述之勢能參數列於表 3.2 中。 3.2.3. 金屬勢能 本研究使用金原子作為模具材料，採用 Morse 勢能函式並選擇 Doyama 等人 43所研究的金原子參數。該形式為較簡單的二體勢能，函式如下所示：

(

)

(

)

(

)

[

exp 2 r r0 2exp r r0

]

D U_Morse = − α − − −α − ( 3.9 ) 其中D為勢能能量參數。上述之勢能參數列於表 3.3 中。 3.2.4. 其他分子間勢能 對於上述的三種材料之間的交互作用力，能否正確的使用適當的勢能及參數，對於本研究有十分重大的影響。聚氧乙烯之所以選擇計算量較大的 EA model 主要是考量到與水之間交互作用的描述，Smith 等人36_{便針對這部分進行研究，並提供了與各種水分子模型的參數，} 其勢能與參數列於表 3.4。同樣的，金原子與高分子溶液之間的交互作用力也十分重要，Dou 等人44_{提供了金與水之間的勢能，如下所示：}

(

)

(

)

(

)

[

1 1

]

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

(

)

(

2

)

0 exp 2 H e H Au D r r U ₋ =γ − β − ( 3.11 ) 該函式由 Spohr45所提出，與 Morse 勢能相似。考量到計算效率，因此研究中會使用到截斷勢能，並使用開關函數(Switch function)46_進

(38)

行修正，讓勢能函數能平緩的衰退並於截斷距離歸零，該函式如下所示： ( ) ( ) ( )

(

) (

)

(

)

⎪ ⎪ ⎩ ⎪⎪ ⎨ ⎧ ≤ < − − + − ≤ = = ₋ − 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 ₂ ₃ 1 ( 3.12 ) 其中S( )r 為開關修正式，r_on為開啟半徑，r_off 為關閉半徑。另一方面，金原子與聚氧乙烯的分子間作用力採用 Lorentz-Berthelot Mixing Rule 獲得，關係式如下所示：

(

A B

)

/2 AB σ σ σ = + ( 3.13 ) B A AB ε ε ε = ( 3.14 ) 使用 Lennard-Jones 12-6 勢能描述，參數列於表 3.3 中。

(39)

表 3.2 Poly (ethylene oxide) 與 1,2-dimethoxyethane 勢能參數表(Smith, Jaffe et al. 1993; Chunxia, Praveen et al. 2006)

(40)

表 3.3 Au 相關勢能及參數(Imafuku, Sasajima et al. 1986)

(41)

3.3. 模擬流程

(42)

3.4. 初始化位置及速度

初始位置的排列並沒有一定的限制，不過通常會選擇最接近物質真實狀態的方式來排列，可使系統較快達到平衡。例如當對象為固體時，會採用接近自然界的晶格排列方式來給予初始位置，若為氣體或液體的系統時，給予隨機的排列則會比較有利。較常見的固體晶格排列有四種堆積方式，分別為簡單立方堆積 (SC)、體心立方堆積 (BCC)、面心立方堆積 (FCC)及六方最密堆積 (HCC)。分子動力學模擬中，噴嘴多以單一金屬建立，所以通常使用 FCC 排列。水的模型則是使用 BC 或是 FCC 排列。另外，高分子的部分只要保持最穩定型態，分子間取得適當距離即可。建立初始位置後，藉由亂數給定每一顆粒子的速度，並將速度分量修正到一個正負範圍之內，使粒子的速度呈現隨機的方向。透過以上的位置及速度的配置，然後經過長時間的平衡，即可使粒子呈現 Maxwell 速度分佈而接近真實的狀態。另一方面，在模擬上通常會將系統之質心速度歸零，使系統不會移動，然而透過這種方式給予速度並不能確保系統總動量為零，因此必需在亂速給定之後將總動量修正，其方程式如下： 1 / N new i i i k k i k m m m N ⎛ ⎞ =_⎜ − _⎟ ⎝

∑

⎠ v v v ( 3.15 ) 其中N為總粒子數，m為質量，v為粒子速度。

(43)

3.5. 系統類型與控制

在分子動力學中，會依照現實中的系統選定需要的系統類型，常見的類型分別為 NVE ensemble、NVT ensemble 及 NPT ensemble。 NVE ensemble 即粒子數 (N)、體積 (V) 及內能 (E) 為恆定值，在奈米噴流的系統中先進行 NVT ensemble 獲得噴流前系統，然後在不控溫的條件下進行噴流模擬。

3.5.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 = > = − < = + ( 3.16 ) 2. 最小鏡像法：為了描述巨觀性質，當粒子間距離超過截斷半徑時，必需將粒子間距修正到最短距離。目的是為了讓中心粒子周圍如充滿粒子的巨觀系統，但實際上粒子數並沒有改變，所以當某一維度使用週期邊界時，可以視該方向有無限延伸的複製系統所組成，若三維度皆使用，則如同在一個巨觀系統進行計算。

(44)

∑

+ = < − = > 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 ( 3.17 ) , 2 x ij ij ij x L when x < − x =x +L 另外還有虛擬牆及實體牆邊界，在奈米噴流模擬中，噴嘴本身就是以真實粒子組成的實體牆，而推動粒子噴出的驅動力通常也是來自真實推板的推擠。 3.5.2. 系統控溫

1. Simple Velocity Scaling :

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

(45)

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 ( 3.20 ) ζ 為摩擦係數，g為系統自由度，Q為溫控質量，τ 為特徵時間。在該方法下，系統周圍如同存在一個恆溫裝置。當系統內的溫度過高時，則將能量透過摩擦係數釋放能量，反之則提高能量，因此可以使系統在控溫過程中維持能量恆定。

3.6. 數值演算法

在分子動力學中最常見的兩個演算法為 Gear’s 五階 predictor-corrector 及 velocity verlet 兩種數值方法，本研究使用 Gear’s 五階 predictor-corrector，因為該方法可以獲得較精準的預測。演算的過程分成三個部分，分別為預測 (Prediction)、計算 (Evaluation) 和修正 (Correction)，流程如下： 1. 預測 (Prediction) 預測的部分使用 Taylor-series 五階式進行演算，在初始化之後，藉由已知的數值預測出下一個時間點的位置、速度及加速度等。方程式以矩陣方式列於下方：

(46)

( 3.21 ) 2. 計算 (Evaluation) 藉由預測獲得的粒子位置配合勢能方程式計算出每顆粒子所受到的力，並將作用力轉換成為加速度。 3. 修正 (Correction) 由於計算值的正確性要比預測所獲得的加速度可靠，因此將加速度修正為較正確的值，並利用兩個加速度的差值配合適合的參數進行誤差修正。 ( ) ( ) ( ) ( )2 i, 2 e i, 2 i, 2 i, 2 -2! t r r r r R = Δ Δ Δ ≡ Δ ( 3.22 ) ( 3.23 ) 式中的a0、a1、a2、a3、a4及a5為 Gear’s 五階 Predictor-Corrector 的常數值，會隨著演算法的階數 (q) 不同而改變，如下表所示： (1) (1) 0 (2) 2 (2) 2 1 2 (3) 3 (3) 3 3 (4) 4 (4) 4 4 5 (5) 5 (5) 5 1 1 2! 2! 1 1 3! 3! 1 1 4! 4! 1 1 5! e 5! r r r t r t r t t _a a r t t a a r t t a r t t α ⎡ ⎤ ⎡ ⎤ ⎢ _Δ ⎥ ⎢ _Δ ⎥ ⎢ ⎥ ⎢ _{⎥ ⎡ ⎤} ⎢ ⎥ ⎢ _{⎥ ⎢ ⎥} Δ Δ ⎢ ⎥ ⎢ _{⎥ ⎢ ⎥} ⎢ ⎥ ⎢ _{⎥ ⎢ ⎥} ⎢ _Δ ⎥ ₌⎢ _Δ ⎥₊_{⎢ ⎥}_Δ ⎢ ⎥ ⎢ _{⎥ ⎢ ⎥} ⎢ ⎥ ⎢ _{⎥ ⎢ ⎥} ⎢ _Δ ⎥ ⎢ _Δ ⎥ _{⎢ ⎥} ⎢ ⎥ ⎢ _{⎥ ⎢ ⎥}_{⎣ ⎦} ⎢ ⎥ ⎢ ⎥ ⎢ Δ ⎥ ⎢ Δ ⎥ ⎢ ⎥ ⎢ ⎥ ⎣ ⎦ ⎣ ⎦ r R r r r 2 3 4 5 2 3 4 (1) ( 2) 2 3 (3) ( 4) 2 (5) 1 1 1 1 1 2! 3! 4! 5! ( ) ₁ ₁ ₁ 0 1 ( ) _2! _3! _4! ( ) 1 1 0 0 1 2! 3! ( ) 1 ( ) ₀ ₀ ₀ ₁ 2! ( ) 0 0 0 0 1 0 0 0 0 0 1 t t t t t r t t r t t t t r t t r t t t t t r t t r t t _t _t r t t t ⎡ _Δ _Δ _Δ _Δ _Δ ⎤ ⎢ ⎥ ⎢ ⎥ + Δ ⎡ _{⎤ ⎢} _⎥ Δ Δ Δ Δ ⎢ _{+ Δ} ⎥ _⎢ _⎥ ⎢ _{⎥ ⎢} _⎥ ⎢ + Δ ⎥ _⎢ _⎥ Δ Δ Δ = ⎢ _{+ Δ} _{⎥ ⎢} _⎥ ⎢ _{⎥ ⎢} _⎥ ⎢ + Δ ⎥ _⎢ _Δ _Δ _⎥ ⎢ _{⎥ ⎢} _⎥ + Δ ⎢ ⎥ ⎣ _{⎦ ⎢} _⎥ Δ ⎢ ⎥ ⎢ ⎥ ⎣ ⎦ (1) (2) (3) (4) (5) ( ) ( ) ( ) ( ) ( ) ( ) t r t r t r t r t r t ⎡ ⎤ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎣ ⎦

(47)

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

表 3.5 三~五階 Gear Predict-Corrector Algorithm 參數

3.7. 加速計算方法

3.7.1. 牛頓第三運動定律 由於分子動力學假設分子運動符合牛頓運動定律，當使用作用力與反作用力時，計算兩粒子間彼此的作用力時僅需要計算一次即可，兩粒子只有受力方向不相同而已，因此可以減少一半的運算量。 3.7.2. 分子勢能截斷法 計算勢能以及粒子間作用力時，根據勢能方程式我們可以知道當粒子間之相對距離夠遠時，所獲得的數值非常小，因此，我們會設定一個截斷半徑 (Truncated potential)，當距離大過該半徑時便可將忽略。換句話說，只要使用該方法，在超出範圍外的粒子可以不必計算，如此一來我們將可節省大量的運算資源。

(48)

( )

r r r r r r r U ij c ij ij ij ij > ≤ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ − ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ = , 0 , 4 6 12 σ σ ε _{( 3.24 )} 除此之外，使用截斷半徑除了節省粒子間交互作用的運算量，更重要的是可以配合 Neighbor List 及 Cell-Linking 法，更加的提升運算效率。然而，使用截斷半徑將引發其他的問題及限制條件，分為以下三點說明： 1. 作用力平移 (shifted-force potential)：當截斷勢能後，方程式將出現不連續點。根據漢米頓力學，分子間作用力為勢能方程式微分並且取負值，因此截斷勢能將會造成誤差，因此提出作用力平移修正不連續點，方程式如下：

( )

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

( )

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

∫

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

(49)

圖 3-2 勢能 LENNARD-JONES 12-6 SHIFTED POTENTIAL 示意圖 2. Long-range correction: 雖然在截斷半徑之外的勢能很小，但是若使用週期邊界描述巨觀性質，則截斷半徑之外應該有非常龐大的粒子數被忽略，僅管半徑之外勢能很小，但是若粒子數非常多時將會造成很大的誤差。Long-range 修正便針對這點提供合理的修正。 ( ) ( )

∫

∞ = C r C _U _r _g _r _r _dr N U ₂ 2πρ ( 3.28 ) 其中ρ為密度，g( )r 為徑向分佈函數。 3.7.3. 鄰近列表法 分子動力學的計算量非常的龐大，僅管已經使用了截斷勢能法，運算量隨著粒子數目的半方成正比關係，因此提出鄰近列表法 (Neighbor-list method) 加速計算。

(50)

鄰近列表法的基本觀念，在系統中的任一分子，其周圍的分子於短時間內變動並不會太大，意即於此段時間內，此分子僅需和周圍鄰近固定的分子計算作用力，可以不用對整個系統的分子作計算；當此方法使用在越大的系統時，加速的效果更為顯著。圖 3-3 鄰近列表法示意圖鄰近列表法的使用方法，先將分子 i 列表半徑(List Distance ,r1) 內的分子都列入清單中，在接下來的時間間隔中，分子 i 只需對清單內的分子作距離計算，持續到下一次的更新列表，而鄰近列表清單只需要數個或數十個時間間隔後再予與更新。其中列表半徑和清單列表的頻率是相互影響的，選擇的決定依據在於，位於列表半徑外的分子能否在下次清單更新前，進入影響分子的截斷距離內，一般而言，當分子的移動速度越快，相對的列表更新頻率增高或者是列表半徑要增加。 3.7.4. Cell-linking method 此方法使運算量僅與分子數目成一次方關係，大量減少運算量，並且可與鄰近列表法做結合，更有效率的節省計算時所花費的時間。此方法的概念是先將模擬的系統切割成若干等大小的方格(Cell)。

(51)

判斷出系統中每個分子所處的方格，當欲計算系統某一個分子受力情形時，只需找尋在分子所屬方格周圍附近，未超過參考距離的數個方格內的分子來判斷即可，不需要對整個系統的分子個別做搜尋判斷，而考慮使用的距離應為截斷距離，但絕大部分此方法都會和鄰近列表法同時使用，則參考距離就採用列表半徑。圖 3-4 CELL-LINKING 法示意圖

3.7.5. 列表勢能法 (tabulated potential method)

本研究中水與高分子溶液的結構特性研究當中必需使用 Ewald sum 計算庫侖力；高分子溶液噴流研究中使用到 Morse 及 Lennard-Jones exp-6 勢能，皆很容易遇到三角函數、exponentiation 及 error function 等函數，若每次都要對這些函數進行計算勢必會大大的耗費許多時間。若將這類函式或勢能和力直接列表，模擬時以查表方式讀取數值，則可以省去很多的計算時間。因此本研究中把 erfc、exp 及三角函數皆列表運算，並分割為數萬等分，來節省模擬的時間。

(52)

3.8. long-range interaction

勢能的運算中，包括 Lennard-Jones 及 Morse potential 皆屬於 short-range interaction，然而當勢能屬於 long-range interaction 時，一般的運算方式並不能正確的描述，本節將討論如何處理 long-range interaction 的運算。在short-range interaction計算塊狀性質，已經普遍的使用 3.5.1 提及的週期性邊界條件，然而為了避免不正確的表面效應，因此允許使用 3.7.2 所提到的截斷勢能法同時也帶來了提高運算效率的好處。然而對於long-range interaction而言，相同的處理方式將會發現到截斷空間之外的貢獻十分的大，不合理的截斷將造成極大的誤差及不自然的分子運動。當勢能允許使用截斷半徑時，都必需可以使用 long-range correlation補償截斷空間之外被忽略的貢獻。假設存在一勢能函式為 d r− 形式，並代入 3.7.2 的long-range correlation的勢能修正式中並積分，當d ≥3時，修正值將為無窮大，意即無法使截斷勢能。對於這類分子間相對距力增加時，勢能之交互作用降低很慢者視為long-range

interaction。之前提到的庫侖力 (charge – charge, −1

r ) 的運算即屬於

long-range interaction。

若將庫侖力與凡得瓦力勢能相互比較，如圖 3-5 所示。在截斷距離處可以看到 Lennard-Jones 幾乎為零，而庫侖力明顯存在相當大的數值，截斷之外的勢能仍會對系統造成相當明顯的影響。

(53)

圖 3-5 比較庫侖力與 LENNARD-JONES 勢能降低至零的程度 3.8.1. Ewald sum 庫侖定律視電荷為點電荷進行運算，然而我們知道上述的系統若直接運算並不實際。在 Ewald 的計算法中，將每一個點電荷周圍附加正負相反，大小相等的電荷分佈，並運用較合宜的高斯分佈來描述。點電荷在覆蓋如同離子大氣一般，遮蔽與鄰近電荷的交互作用，經推導後為 short-range interaction 的 error function 形式。

然而在額外添加電荷分佈後，必需加入抵消(canceling)的電荷分佈項來進行補償。而該抵消項經過傅利葉轉換至虛數空間(reciprocal space)累加後，再將總值轉換回到真實空間(real space)求得其貢獻值。另外，由於虛部空間項的累加中將粒子自身的作用也一併累加進去，因此必需額外再減去一項 self-term，最後可將所有方程式整理如下

(54)

(

)

(

)

∑∑

∑

= = = = ≠ ∞ = ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ + ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ − ⎟⎟ ⎟ ⎟ ⎟ ⎟ ⎠ ⎞ ⎜⎜ ⎜ ⎜ ⎜ ⎜ ⎝ ⎛ • ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ − ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ + + + = N i N j N i N i i i i ij j i ij ij j i coul q L q k q q L erfc q q U 1 1 1 2 1 3 2 2 / 1 0 2 2 3 0 ' 3 2 cos 4 exp 1 2 1 r r k n r n r k n π π κ κ π κ ( 3.29 ) 其中κ值愈大將使真實空間項的 error function 降至零的速度愈快，也會使高斯電荷分佈的波峰提高，這點將會影響到虛部空間應該取多少個 k-space。一般會選擇讓真實空間能在n=0即收斂至接近零，因此運算上將保持和一般的週期鏡像的計算方法相同，而在虛部空間取適當的 k-space 進行累加。實際運算上，一般取κ =5/L，真實空間的截斷半徑取L/2，虛部空間的 k-space 約為 100~200 進行累加。注意到式中的第四項，該項考慮球形系統之外的介質對系統的影響。當環境介質為良好的導體時無需添加該項，若為真空環境則必需要進行修正。

3.9. 性質計算

3.9.1. 溫度 以微觀的角度來看，溫度 (Temperature) 來自於原子的運動，分子動力學可以輕易的追踨每顆粒子的位置及動量，籍此計算出系統的溫度。根據氣體動力論可以得到下列方程式：

∑

= ⋅ = = N i i i B i B k N k m T T k E 1 3 2 3 v v ( 3.30 )

(55)

其中Ek為動能，N為總粒子數，kB為波茲曼常數，T為溫度。 3.9.2. 壓力 壓力 (Pressure) ： ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ + = + =

∑

∑∑

< = i j ij ij N i i i fx mx m V P P P v r DF 1 2 3 1 ( 3.31 ) 其中Pmx為粒子通過空間平面之動量變化量造成的壓力。Pfx為兩粒子間交互作用施加在平面上的壓力，Fij為 i 與 j 粒子間作用力。壓力的 long-range correction 列於下： ( ) ( )

∫

∞ − = c r B B LR _g _r_r _dr dr r dU r T k T k P 2 3 2πρ ρ ( 3.32 ) 其中ρ為密度，kB為波茲曼常數，T為溫度，U( )r 為勢能函數，g( )r 為徑向分佈函數。 3.9.3. 徑向分佈函數 在分子動力學模擬中無法在模擬前就確定物理的相態，即便是參考了實際的文獻資料，也不能斷定在模擬中會有相同的結果，因此 RDF 通常會被拿來做決定結構的方法。

( )

(

)

(

r r

)

V N M r r N r g M k k Δ ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ Δ =

∑

= , 2 1 , 1 ρ ( 3.33 ) (r r) N_k ,Δ 為每個 time-step 在距離r~r+Δr球殼中出現的粒子數，N 為總粒子數，ρ為系統平均密度，V(r,Δr)為距離r~r+Δr球殼體積，M

(56)

為統計 M 步 time-steps。其物理意義為一個粒子質心出現在另一個粒

子質心r~r+Δr球殼空間的機率，在奈米噴流中可以藉由徑向分佈函

(57)

4. 模擬系統

4.1. 水與聚氧乙烯水溶液塊材系統

於研究的第一階段先建立高分子溶液以作為噴流模擬的材料，在選定合適的勢能之後，將分為水溶液與聚氧乙烯水溶液兩系統進行研究。水溶液的勢能主要是根據Yujie等人於 2006 年發表之文獻 22_內容來建立系統，系統為NVT ensemble，使用 216 顆水分子，並固定體積保持密度為 3 / 1g cm ，平衡與取樣的模擬採用Nose-Hoover控制系統溫度，使溫度維持在 298.16 K。time-step設為 0.4 fs，模擬總共執行 2 ns (5000000 time-step)，1 ns使系統達到平衡，再執行 1 ns 進行水溶液的驗證。該系統模擬條件詳列於表 4.1 中。於聚氧乙烯水溶液的部份，根據Smith等人 33 ~36 , 48 _{一連串的研} 究，可以瞭解通常於分子動力學模擬中都會以 1, 2-二甲氧基乙烷 (DME)來代替聚氧乙烯，主要是因為 1, 2-二甲氧基乙烷可以說是聚合度為 1 的聚氧乙烯，且所有的勢能參數皆與聚氧乙烯相同。因此針對 1, 2-二甲氧基乙烷水溶液進行研究及討論，其性質即可代表聚氧乙烯與水之間的性質。包括聚氧乙烯與水之間的作用力，都是先透過 1, 2-二甲氧基乙烷的研究建立參數，進而延伸至高分子聚氧乙烯的相關研究。模擬系統與水溶液同樣為NVT ensemble，使用兩種不同的溶液組成，分別為XDME =0.18,0.72莫耳分率 (分別為 46 條DME加 200 顆水分子及 77 條DME加上 30 顆水分子)。固定體積保持密度為 3 / 122 . 1 g cm ，並控制溫度保持在 318 K。模擬總共執行 0.8 ns (2000000 time-step)，

(58)

0.1 ns使系統達到平衡，再執行 0.7 ns 進行取樣計算。該系統模擬條件詳列於表 4.2 中。

圖 4-1 水分子平衡示意圖

(59)

Model SPC/Fw Ensemble NVT Molecular numbers 216 H2O Density 1 g/cm3 (L is 18.6 Angstrom) Temperature 298.16 K Thermostat Nose-Hoover

Algorithm Fifth-order Gear’s predictor-corrector Long-range interaction Ewald sum

reciprocal space of k 6

Real space of cut-off L/2 (L is system of boundary length) Short-range of cutt-off 9 Angstrom

Total simulation time 2 ns (5,000,000 time-steps) Equilibrium time 1 ns (2,500,000 time-steps) Sampling time 1 ns (2,500,000 time-steps) Time-step 0.4 fs

表 4.1 研究塊材水溶液性質之系統操作條件 Model Explicit atoms (by Smith et al.) Ensemble NVT

Molecular numbers 30 DME & 200 H2O XDME =0.18 77 DME & 30 H2O XDME =0.72

Density 1.122 g/cm3 (L is 23.5 Angstrom) Temperature 318 K

Thermostat Nose-Hoover

Algorithm Fifth-order Gear’s predictor-corrector Long-range interaction Ewald sum

reciprocal space of k 6

Real space of cut-off L/2 (L is system of boundary length) Short-range of cutt-off 9 Angstrom

Total simulation time 0.8 ns (2,000,000 time-steps) Equilibrium time 0.1 ns (250,000 time-steps) Sampling time 0.7 ns (1,750,000 time-steps) Time-step 0.4 fs

表 4.2 研究塊材 DME 水溶液性質之系統操作條件

4.2. 高分子溶液噴流系統

(60)

圓柱形狀的噴嘴系統。噴嘴大小為 9.168 nm、內部直徑為 6 nm ，出口直徑為 3 nm，進出口截面積比為四比一，噴嘴收縮角度為 90°。詳細的幾何結構如圖 4.2 中介紹。圖 4-3 嘴模型幾何示意圖噴流所使用的高分子溶液製備系統上則使用 6240 顆水分子及 6 條聚合度 12 的聚氧乙烯分子。在進行噴流模擬之前，必需將材料控制到指定的初始條件然後放入噴嘴之中。首先將材料適當的排列在適當大小的空間中，透過一作用力將材料壓縮到指定體積，然後將材料置入噴嘴中，並進行 0.2ns 使材料與噴嘴達到平衡狀態。在平衡過程中，噴嘴內部高分子溶液及內緣金屬層控制在指定溫度，由於噴嘴表面在噴流過程中將產生潤溼(wetting)的現象，因此金屬溫度需略高於溶液的溫度，可有效降低潤溼現象的情形。材料部分使用的是 Nose-Hoover 進行控溫，噴嘴部分考量到執行噴流過程為非平衡態，因此使用 Simple velocity scaling 控溫。開始執行噴流時，將不再控制材料的溫度於真空下進行模擬。

(61)

5. 結果與討論

5.1.SPC/Fw 勢能模型之分子結構 首先進行水溶液系統進行的模擬並與文獻22_{的結果比較。根據文} 獻22的研究，提出調整平衡鍵長及鍵角參數可以改變SPC model某些特定的統計性質，同時些許影響水分子結構，本研究將透過徑向分佈函數與文獻值進行比對。模擬系統使用與文獻相同 216 顆水分子，於 298.16 K，密度為 3 / 1g cm，平衡鍵長與平衡鍵角分別為 1.012Å及 113.24 °，詳細系統條件列於表 4.1 中，研究中統計氧-氫原子、氧-氧原子的徑向分佈函數。其徑向分佈函數模擬結果如圖 5-1 所示，與文獻22 結果幾乎完全相同，如圖 5-2。可以發現該分子模型的模擬結果與實驗值相當接近，僅僅只有在第一水合層的部分有少許高估的情形。

(62)

圖 5-1 水分子 216 顆，溫度 298.16 K 密度為 3

/

1g m ，SPC/FW MODEL 徑向分佈函數之模擬結果圖。

(63)

(64)

5.2.高分子溶液結構性質

成功的驗證水分子之後，更進一步研究高分子溶液材料。如同前面所提到的，1, 2-二甲氧基乙烷(DME)可以代替聚氧乙烯與水溶液之間的作用力，因此針對DME水溶液進行研究及討論。研究中分別使

用 45 條DME與 200 顆水分子混合(XDME =0.18)及 77 條DME及 30 顆水

分子混合(X_DME =0.72) 兩種系統，於 318 K密度為 3 / 122 . 1 g m 條件下進行模擬，詳細系統條件列於表 4.2 中。研究中統計氧(DME)-氫(水) 原子、氧(DME)-氧(水)原子的徑向分佈函數及不同莫耳分率下的局部密度比值，並與文獻 36_{的研究結果進行比較。由圖 5-3 不同DME莫耳} 分率下的徑向分佈函數及圖 5-5 不同DME莫耳分率下的密度分佈圖與文獻 36_{幾乎有相同的結果 (圖 5-4、圖 5-6)。由此可得知程式在模} 擬高分子溶液部份為正確的。

(65)

圖 5-3 DME 水溶液於 318 K，DME 莫耳分率(XDME)為 0.18。分別為

ODME-OWATER及 ODME-HWATER的徑向分佈函數。其中 ODME-HWATER的部分補償

2 的數值。

圖 5-4 文獻36_{中 DME 水溶液 O}

DME-OWATER及 ODME-HWATER的徑向分佈