臺灣博碩士論文加值系統

摘要

本論文利用分子動力學的方法(Molecular Dynamics，簡稱MD)模擬水分子在白金觸媒層表面的接觸反應情形。在微尺度之下，液體和固體之間的接觸現象，主要是因為分子之間的作用力而造成，我們利用分子動力學的方法來模擬微燃料電池中水分子和白金團簇之間的相互反應情形，觀察接觸角的變化情形與水分子在白金觸媒層上的擴散情形。

模擬系統的邊界條件為週期性邊界條件、鏡面反射(Mirror Reflection)邊界條件、而白金觸媒層邊界利用Lennard Jones Pt-Pt的勢位能處理，模擬的水分子是利用已知的SPC/E model，而水分子對於白金團簇的勢位能關係是採用E.Spohr根據 所計算出來的勢位能方程式來進行我們的計算模擬，我們初步關心的問題是水對於白金表面的濕潤情形，因此我們模擬水分子液滴在白金表面的反應情形，計算接觸角擴散速率和時間之間的相互關係以及接觸角的變化情形。

Abstract

This thesis utilizes the molecular dynamics method simulation water molecules in contact with platinum catalyst surface. Under microscale, the contact phenomenon between the liquid and solid, mainly cause of effort between the molecules. We use the molecular dynamics method to simulate the water droplet in contact with platinum catalyst surface in the micro fuel cell. We also observe the change situation of contact angle and water spreading velocity on the platinum catalyst surface.

The simulation system have periodic boundary conditions and mirror reflection boundary condition, we use Lennard Jones potential for platinum catalyst molecules. The water molecules of simulation utilizes known SPC/E model, and we use S-H potential between water molecules and platinum molecules, it is calculated by E.Spohr and Huckel. We care about tentatively is the moist situation of water on the platinum surface , so we simulate water droplet on the platinum surface, and calculate the change situations of the reaction of contact angle and water spreading velocity on the platinum catalyst surface.

第一章 緒論

1.1引言

隨著科技日新月異的進步，每年全球汽機車量逐年成長，汽機車的動力來源為內燃機，傳統內燃機一直存在著空氣污染的問題，內燃機的燃料為石油，石油是地球有限資源，在人類過渡使用的情形下，全球石油的蘊藏量終將有使用完的一天。有鑑於此，科學家不斷積極的在研發低污染、高效能的新能源技術，燃料電池是新一代能源技術的里程碑，在低溫型燃料電池中，直接甲醇燃料電池(DMFC)、質子交換膜燃料電池(PEMFC)，電化學反應的核心部份是膜電極組 (membrane electrode assembly，簡稱MEA)，其結構可細分為氣體擴散層和觸媒層，氣體擴散層的主要功能是讓反應氣體到達觸媒層，同時防止電解液或水的逆滲透，氣體擴散電極大部分為多孔性的電極，一般為碳紙或碳布組成，觸媒和液體接觸的情形將直接影響其效能，模擬水分子和甲醇分子在觸媒層上的反應情形，此部分為本論文切入的重點。

當尺寸到達微尺度之下(microscale and nanotechnology)，傳統巨觀角度去分析已不再適用。在電腦晶片技術的進步下，分子動力學的方法越來越受到重視，在直接甲醇燃料電池(DMFC)中，水、甲醇和觸媒層的反應，不管是在陰極或者是陽極的作用情形都將影響整個燃料電池的整體效能，在本論文中，初步先架構水分子和觸媒的反應，甲醇的部分由於勢位能參數搜尋困難的問題，暫時沒加入計算，理論上液體和固體之間的接觸現象主要是因為分子之間的作用力而造成，這裡利用分子動力學的方法模擬水分子在白金表面的反應情形。分子動力學法(molecular dynamics，簡稱MD)之基本原理是假設系統中個別一分子(或原子)的運動行為遵守牛頓第二運動定律，總合能量則由分子的位能(potential Energy)及動能(kinetic Energy)來分配。分子與分子之間的作用力，取決於一位勢位能函數，此勢位能產生於分子間的凡得瓦力作用、電荷靜電力、磁偶矩、與分子內部勢位能所造成，最常用的勢位能為Lennard-Jones potential (12,6)，對於水分子，則必須找出其適用的勢位能函數(SPC/E)，方能得到準確的預測結果。

1.2 文獻回顧

在微燃料電池中，液體和固體之間的接觸現象主要是因為水分子和白金分子之間的作用力而造成，關於動力問題的計算有稱為分子動力的計算方法，我們利用分子動力學的方法來計算分子運動行為，Haile一書中提供完整的分子動力學模擬架構[1]，M.P.Allen、D.C.RAPAPORT和R.J.Sadus所著的書[2][3][4]也對液體分子的模擬方法有詳細精闢的講解和說明，Shiegeo Maruyama也利用分子動力學的方法研究微小尺寸的熱傳現象[5]。漢米頓（Hamiltonian）理論讓物理學家從與牛頓力學不同的角度(phase-space)來描述物理現象，但這些漂亮及有用的數學在使用上有其基本上的限制，對於實際的物理問題，往往很困難得到解析解答，隨著電腦的發明及電腦功能的快速成長，可以利用電腦數值的計算得到答案。這個方法先算每一個分子所受的力，然後以數值方法解牛頓的運動方程式以預測原子經一小段時間( t)後估計位置及速度，將分子的位置及速度換成新的以後，再計算所受的力，再預測新的位置及速度，這樣就可一步接一步的算出分子的位置及速度與時間的關係，在從位置與時間的關係，可知分子的運動路徑，利用此種方法研究水分子如何在白金表面的反應情形。

1987年Berendsen[6]等學者，修改之前發表的水分子勢位能函數方程式，水分子是屬於極性分子，由於前面所推導的函數計算分子之間極性效應和實驗值誤差較大，所以使得水分子不能產生凝結狀態的性質，因此水分子的勢位能函數改良導入極性效應和偶極矩，最後形成SPC/E(Simple Point Charge Extended)的模式。

1989年E. Spohr[7]建立水分子對於白金團簇的勢位能函數，其勢位能函數模式的基礎是利用 計算水分子和白金團簇之間的相互反應關係，而得到兩分子相互之間勢位能函數的解析解。

1992年J. Nieminen[8]首先利用分子動力學的方法模擬液滴降在一均質的固體基質表面上的反應情形，其結果觀測到液滴在固體表面上的反應為開始時以一定的速度慢慢的在固體表面擴散開來，然後慢慢的降低擴散的速度，整體隨時間在固體表面擴散的情形接近線性函數的關係，而J.Yang[9]也用分子動力學的方法模擬Ar液滴在固體表面上的反應情形，得到趨勢為 。

Matsumoto等[10]在1995年利用分子動力學的方法，勢位能是用Lennard-Jones勢位能函數，模擬數百個氬液滴在平板固體表面的凝結行為，其中固體表面結構是單層面心力方結構的分子，模擬的結果發現親合性是氬液滴和固體表面凝結形狀的主要因素，液滴會逐漸在固體表面散開，這是以微觀的角度來看液滴的反應情形；如果是以巨觀的機制來看，接觸角的形成可以被解釋成表面能量之間的平衡的反應。1997年，Matsumoto等[11]繼續利用分子動力學的方法，模擬氬液滴/氬氣體和固體分子接觸時的熱傳現象，這時候固體表面是用三層面心立方結構的固體分子，結果發現微觀的模擬方式可發現接觸角的餘弦函數是積分反應勢能深度的線性函數；以巨觀的觀點來看，楊式方程式也是餘弦函數的關係式，對於不同溫度的固體表面模擬，發現其接觸角在平衡狀態時幾乎相同，液滴在凝結或蒸發的溫度分佈、速度分佈和準平衡狀態所量測的熱通量相同，證明利用分子動力學的方法來模擬分子之間的反應準確性高。

S.Maruyama和T.Kimura利用分子動力學的方法作了一系列的研究，在2000年[12]模擬Lennard-Jones分子的蒸汽氣泡在固體表面的成核(nucleation)情形，固體表面的處理方式是用三層結構的幻想分子概念處理，目的是為了維持溫度不變。接著在2001年[13]，模擬水分子在白金團簇上的反應情形，其中水分子利用以知的SPC/E模式，白金表面是用一層的2900顆白金分子，結構為面心結構<111>排列，而水和白金之間的勢位能用E. Spohr所導出來的S-H勢位能，結果發現當水分子越來越多時，水分子在白金表面的擴散情形越來越趨近於 或是 的關係式，其中R為水在白金表面擴散面積的半徑，其中A為水在白金表面擴散面積，而利用L-J勢位能模擬的預測結果大約是趨近於 ，至於實驗所得到的結果是趨近於 。2003年[14]利用相同的理論基礎，將水分子和白金團簇的反應加入由Zhu和Philpott所推導的Z-P勢位能關係，比較S-H和Z-P兩種勢位能的強度以及在不同面心立方結構<111><100><110>表面的反應情形，研究結果得到水分子在白金表面擴散情形一開始趨近於 的關係，然後變成 的關係式，此外水分子的反應結構在Z-P勢位能比較強，而在<100>結構上的接觸角是最大，在<110>結構上的接觸角是最小。

1.3研究的目的和方法

在微觀的研究中，傳統巨觀角度去分析已不再適用，取而代之的是分子動力學的方法，在微燃料電池中，液體和固體之間的接觸現象主要是因為分子之間的作用力而造成，因此我們要利用分子動力學的方法，研究水分子和白金團簇之間的相互反應情形，觀察接觸角的變化情形。

本篇論文利用分子動力學的理論方法建構系統進行模擬，我們關心的問題是水、甲醇對於白金觸媒表面的濕潤情形，因此我們初步先模擬水分子液滴在白金表面的平衡狀態，計算接觸角(contact angle)與接觸面積和時間之間的相互關係。模擬的系統為四個邊界(前後左右)為週期性邊界條件，上面邊界為鏡面反射(mirror reflection)的方式處理，底面白金固體邊界利用Lennard Jones Pt-Pt的勢位能來處理，假設白金為固定的分子，模擬的水分子是利用已知的SPC/E model[6]，而水分子對於白金團簇的勢位能關係是採用E. Spohr[7]根據 所計算出來的勢位能方程式來進行我們的模擬。

1. 4論文架構

本篇論文分為四個章節，略述於下：

第一章為序論的部分，簡略說明如何建立起我們所需要模擬的基本架構，說明我們為何要利用分子動力學的方法來模擬水分子和白金團簇之間的反應，並探討近年來的各種相關研究方法。

第二章為分子動力學的基本架構說明，介紹基本的運動方程式，以及我們所使用的相關勢位能函數，討論邊界條件的處理方式，並說明我們計算勢位能函數的方法，最後說明接觸角的計算方式。

第三章為分子動力學模擬的基本流程圖介紹，這章討論我們所建立模式的基本架構，簡略說明我們所使用程式的基本運算流程和計算的方法，然後介紹單位的無因次化，簡單說明我們所使用的預測分子碰撞軌跡的數值方法以及加快電腦運算的方法，最後為轉動所使用的理論簡述。

第四章為結果的討論部分，這章整理我們所模擬得到的結果，討論接觸角的產生情形以及水分子擴散附著於白金觸媒的速度和時間的關係，分析接觸角的物理變化現象。

第五章為結論和建議的部分，這張我們整理之前所模擬得到的結果，整理出幾點結論，並對未來研究的方向提供建議。



第二章 分子動力學理論

2.1理論

分子動力學所涵蓋的範圍很廣，以模擬而言，一般是以適當的理論方法來描述系統中的分子，預測分子的碰撞情形或解釋分子系統的性質(如熱力學性質)。隨著電腦科技的進步，近年來的分子模擬常利用電腦的高速運算能力，模擬較大的分子系統或用較嚴謹的理論方法以求得較精確的預測結果，而藉由電腦將分子模擬的情形和各種性質轉化成視覺上易瞭解的圖形也是常用的方法。

分子動力學理論之基本原理為假設系統中個別一分子或原子的運動遵守牛頓運動定律(Newton’s second law)方程式，分子之運動加速度由其分子質量及分子所受外力之總合決定，而分子的受力，主要來自分子與分子之間以及分子與壁面固體分子之間的相互作用力，而其總合能量則由分子的位能及動能(包含平移、旋轉、及振動等動能)其方程式為(2.1)，分子與分子之間的作用力，來自一位勢能函數，此位勢能係由分子間的凡得瓦力作用、電荷靜電力、磁偶矩、與分子內部勢能造成，當二分子距離甚近時，此勢能對分子距離(rij)的導數將得出一排斥力，而在某適當距離時，則表現出吸引力(圖2-1)。利用MD理論計算分子的電荷分佈、分子內及分子間的力場，提供分子以及更大粒子的多質點動力計算之基礎，此外，由分子或粒子動力計算並透過統計力學的方法，可計算取得物質之熱力性質(分子自我組織化、相圖、熱傳、擴散率、力學性質等)和相介面性質(表面張力、接觸角、邊界熱阻抗、熔融與結晶等)，對於單原子分子(惰性分子)，最常用的位勢能為Lennard-Jones potential[1]，對於其他種類的分子，則必須找出其適用的位勢能函數，才能得到準確的模擬結果，對於水分子，我們採用已知的SPC/E model[6]。
(2.1)
其中
Etotal total energy
Ek kinetic energy
U internal energy

2.2 分子動力學基本運動方程式

分子動力學的主要目的是利用運動方程式來計算分子在空間中的運動行為，計算出分子運動的軌跡，在描述分子的運動行為，可以用兩種方程式來說明，一種是牛頓第二運動方程式，另一種為漢米頓運動方程式，下麵分別詳細介紹兩種運動方程式。

2.2.1 牛頓力學(Newtonian Dynamics)

描述分子運動行為的牛頓第二運動方程式為(2.2)，牛頓第二運動方程式是屬於二階常微分方程式，以三維空間來說可以表示為(2.3~2.5)，如果分子i 不受到外力影響，方程式(2.2)可以簡化成vi=常數，分子的運動為靜者恆靜、動者恆動，也就是牛頓第一運動定律，根據(2.2)，假設孤立系統不受任何外力，則系統總力等於零，因此系統中的分子i和j的力是必須平衡，則可推出牛頓第三運動定律(2.6)作用力等於反作用力。

(2.2)
(2.3) (2.4) (2.5)
(2.6)

其中
i,j atom i and j
Fi force on atom
m mass of molecule
ri postion of atom I relative to a space-fixed origin
Fix the force of x direction
Fiy the force of x direction
Fiz the force of x direction
vi velocity of atom I
t time

2.2.2 漢米頓力學(Hamiltonian Dynamics)

Hamiltonian Dynamics的原理為假設系統中有N個分子，則系統的狀態由分子的位置(rN)和動量(pN)來表示(2.7)，我們定義系統的總能量為Etotal(2.8)，根據方程式(2.7~2.8)，我們可以推導出方程式(2.9)和(2.10)，將之帶入我們利用漢米頓(Hamiltonian Dynamic)運動方程式(2.11)，求得分子運動變化的軌跡。
(2.7)
(2.8)
(2.9)
(2.10)
(2.11)
其中
r position
p momentum
N spherical molecules
2.3 勢位能函數(Potential Function)

勢位能由分子間的凡得瓦力作用、電荷靜電力、磁偶矩、與分子內部勢能造成，其方程式為(2.12)，相鄰兩分子之間距離(rij)，在某距離之內表現出排斥力(repulsive force)，而在某適當距離時，則表現出吸引力(attractive force)。對於不同的分子需選擇不同的勢位能函數，才能使分子動力學的模擬結果更準確，下麵將針對幾個我們模擬中使用到的勢位能做介紹Lennard-Jones potential、SPC/E和S-H potential function。
(2.12)
where
effective pair potential
rij the distance between molecules i and j

2.3.1 Lennard-Jones potential

為了使模擬更準確，不同的分子必須選擇適合的勢位能函數，目前最廣泛使用的勢位能函數為Lennard-Jones (12,6)，其方程式表示如下(2.13)，對於不同的惰性分子，例如：He、Ar、Kr和Xe，有不同的 、 ，分子之間的potential關係如圖2.1，當兩分子距離為 時，勢位能為零，兩分子之間最小的勢位能在(2.13)的一次導函數等於零時有最小之值，也就是當兩分子之間的距離等於 ，將此值代入(2.13)可得到-ε，也就是分子之間最小的勢位能，在這裡我們利用此勢位能方程式來計算氧原子與氧原子之間以及和白金分子與白金分子之間的勢位能，其參數列於(表一)。

(2.13)
其中
energy at minimum of Lennard-Jones pair potential(J)
length scale (nm)

2.3.2 水分子勢位能(Effective pair potential for water)

水分子和原子不同的地方在於真實的水分子是極性分子，會有移動、轉動、和震動等現象，水分子的結構圖如圖2-2，其勢位能方程式為(2.14)，其中前面的式子表示O-O之間的勢位能函數，後面的式子表是為氫原子與氧原子之間的庫倫力，假設有兩個水分子,則其之間共有九組相互影響的庫倫力，水分子模式(SPC/E)參數列於(表一)，一個水分子的結構式是呈V字型，H-O-H夾角109.47度(2.15)，其中O-H是由氧原子和氫原子共用電子所形成的共價鍵(covalent bond)，但是氧原子對共用電子的吸引力較大，所以共用電子(shared electrons)靠近氧原子這端的時間較氫原子長，使得氧原子稍微帶過剩的負電，而氫原子稍微帶正電，由於這種不均等的電荷分配(unequal charge distribution)，所以水是極性分子。1980年代Stillinger及Rahman 所提出的ST2模式是水分子位元勢能的標準式，同時，各種改良的位勢能紛紛被提出，其中最常被使用者有Jorgensen等所提出的TIP4P模式與Carravetta及Clementi所提出之CC模式，目前較常使用的模式為SPC(Simple Point Charge)和SPC/E(Extended SPC )[6]，這些模式主要反應水分子的極性，以及O-O及O-H原子間的作用力及電荷力，對於模擬水的行為有較準確的模擬結果，在計算水分子運動的部分我們只考慮水分子的移動和轉動，震動部分忽略不計。

(2.14)
(2.15)
其中
oo pair potential for water
oo the diameter of the water molecule
R12 the distance between oxygen and oxygen
o 8.852x10-12 the permittivity of free space

2.3.3 表面和水分子之間的勢位能(S-H potential)

在水分子和白金表面之間的勢位能函數是由 所推導出來，1988年時，Spohr和Heinzinger所建立的模式[7]，其方程式如下:






其中
a1=1.8942x10-16J, b1=11.004nm-1
a2=1.8863x10-16J, b2=10.966nm-1
a3=10-13J , b3=53.568nm-1
a4=1.742x10-19J , b4=12.777nm-1
c=11.004nm-1
the length of the projection of the distance vector onto the surface plane

2.4邊界條件(Boundary Conditions)

在MD的模擬中，一般而言系統包含數百到數千的分子，控制系統受到表面效應的影響，也就是分子和系統壁面的相互影響，因此模擬系統壁面和液體之間，系統因該給予液體在接觸壁面時的反應資訊，也就是給予適當的邊界條件。在本論文的模擬中，我們用到週期性邊界條件(periodic boundary condition)、鏡面邊界條件(mirror boundary condition)、固體表面邊界條件(solid surface boundary condition)，下面分別介紹這三種邊界條件。


2.4.1週期性邊界(Periodic Boundary Condition)

邊界條件的處理，在平衡系統的模擬中，大多以週期性邊界的觀念處理，以下簡稱PBC，使用週期性邊界條件來移除表面效應，可以避免邊界條件的設定而改變取樣的均質性，即或系統邊界中存在有固體壁，該固體壁皆以鏡面反射來處理，或以絕熱壁視之，此二方式皆無考慮固體壁對介質分子的熱傳。

使用PBC來模擬系統，使系統保持固定數目的原子，其原理為:假設體積V，系統內包含N個原子，我們假設體積V只是系統中的一小部分，我們稱它為primary cell。我們先考慮二維的情形(圖2-3)，中間的部分為primary cell而周圍有八個重複cell，我們稱為image cell，其大小形狀皆相同，而每個image cell中也包含N個原子，如同primary cell一樣，整個結構系統我們稱它為primary frame，當原子由 移到 ，我們可以表示為 ，當原子從primary cell移動到左下角，可以表示為 ，其中r為向量。以三維的空間來看，每個系統就有27個cell，當原子移動之後向量可表示為(2.16)，當原子在primary frame，我們將他位置表示為(2.17)

(2.16)
其中
L the scalar length of one edge of a cubic cell

(2.17)

其中
i atom i

在MD模擬之中，我們需要的只有primary cell中的原子，因此在三維系統中，簡單轉換立方晶體中的原子位置可表示向量為(2.18)

(2.18)

將(2.16)(2.17)代入(2.18)可以得到(2.19)
(2.19)

因此，在計算cell中原子的位置時，我們只要知道原子在primary cell中的位置再加上向量

週期性邊界條件可以除去不必要的邊界影響，而且PBC對系統的勢位能以及作用力的影響也可以忽略不計，PBC是分子動力模擬中常用的邊界條件。

2.4.2固體表面邊界條件(Solid Surface Boundary Condition)

在白金分子的部分，我們係假設白金分子團簇為固定分子(body fixed)，也就是底面的固體邊界為固定之邊界，對於白金分子的勢位能，我們利用Lennard Jones(12,6) potential 來計算其分子與分子之間的勢位能函，再將所計算的值加入至系統之總能。

2.4.3鏡面邊界條件

因為系統底部是白金團簇所組成的觸媒層，因此再處理系統上端的部分不能使用週期性邊界條件來處理，在這裡我們使用鏡面邊界條件來處理系統上端部分的邊界問題，鏡面反射的原理是分子的運動路徑乃經由Gear’s predictor correction 預測分子下一時間的位置，然後利用下一時間的位置，反算其經過反射之後的位置，經過鏡面反射的分子，其速度和原來預測之分子的速度大小相等，方向相反。

2.5截斷勢位能(Truncated Potnetial)和長範圍修正勢位能(Long-Range Correction)

在擁有N原子的系統中，根據方程式(2.12)，會有(1/2)N(N-1)的原子相互反應，因此在模擬中，相互反應次數會隨著原子的增加而成平方的增加相互反應。更進一步說，如果我們允許相互反應的原子的範圍增加，r到r+ r，則相互反應的樣品總數為(2.20)內當相互反應因此當原子數量越多時，相對的會大量提高電腦運算的時間，因此我們要對我們的勢位能計算方式截斷然後再修正，一般而言我們取rc =2.5 ，在水分子的模擬中我們取rc =4.73 。

(2.20)
其中
the number of atoms in spherical shell of radius and thickness r
the volume of shell

利用截斷勢位能可以減少大量的電腦運算時間，但是被截斷的勢能會對於性質上的計算有某種程度上的影響，並且造成分子之間勢位能的階梯差，因此模擬的結果只能代表性質的一部份，例如內能和壓力都有其相關的勢位能，因此在計算時必須再加入長範圍的修正(Long-Range Correction)，氧原子與氧原子之間的LRC利用(2.24)，而庫倫力的部分我們採用Reaction field method的方法來修正。在MD模擬中，性質的計算基礎為統計熱力學，而統計熱力學中內能的計算為勢位能乘上徑向分佈函數g(r) (radial distribution function)積分而成(2.21)，以截斷長度為為界，則(2.21)可以分成兩部分(2.22)，當 ，假設g(r)=1，於是我們可以把 寫成(2.24)，將L-J potential帶入可以得到(2.25)。

(2.21)
(2.22)
(2.23)
(2.24)
(2.25)
其中
time average of internal energy
E energy
LR Long range correction

2.5.1反應場Reaction Field Method

在計算長距離修正的部分，氧原子與氧原子的長距離修正前面已經提過，而庫倫力的長距離修正有兩個方法Ewald Summation 和 Reaction Field Method，由於Ewald Summation所需要消耗的電腦計算時間過於龐大，我們這裡採用反應場的方法來取代長距離的修正問題，反應場的方程式為(2.26)，反應場的原理為計算有限半徑rc之內形成的spherical cavity，計算其中的電子反應，而在cavity之外的部分則假設為絕緣狀態，我們只計算局部範圍內的電場反應，另外長距離修正庫倫力還需利用weighting factor f(rij)來修正值以期達到更準確的計算結果(2.27)，


(2.26)
(2.27)

Where
1 dielectric constant outside the cavity (7.0)
rt =0.95rc

2.6 最小影像準則(Minimum Image Criterion)

在N個分子的系統中，分子之間的力計算可以由內能對距離的微分得到(2.28)，根據週期性邊界條件的理論，我們可以將(2.28)改寫為(2.29)，因此我們只需計算primary cell中的分子相互作用力，三維系統中可以表示為(2.30)：

(2.28)
(2.29)
(2.30)

利用上式我們可以計算出整個系統原子的作用力，然而我們計算的力為primary cell中的總力，所以如果勢位能在截斷在rc (1/2)，我們必須算和27個cell之間可能性的最小影響距離，所以我們要定義最小影響距離(minimum image criterion)來計算力。以一維系統來說，我們所關心的為αx的三個值，最小距離為|xij-αxL|，有三種可能發生的情形如圖2-5，分別為A、B、C三種可能性。

在情況A中，( )因為primary cell中的i原子和image cell中的j(-1)原子距離大於1/2L(rc)，我們只需討論primary cell中的i原子和j原子的相互作用力即可，因此可以利用(2.30)求得其作用力。

在情況B中，( )， <(-1/2)L(rc)，primary cell中的i原子和j原子不反應，但是i原子和image cell(-1)中的j(-1)原子反應，其距離為 。

在情況C中，( )， >1/2L(rc)，primary cell中的i原子和j原子不反應，但是i原子和image cell(+1)中的j(+1)原子反應，其距離為。


2.7 接觸角(Contact Angle)

固體、液體、和氣體的反應現象在熱傳的相變化裡扮演一個重要的角色，當液體潤濕固體表面時，原本氣-固的介面被液-固的介面所取代，而氣-固與液-固之介面張力的差，稱之為濕潤張力。當氣-固的介面張力大於液-固的介面張力時，也就是固體和液體間的吸引力大於固體和氣體間的吸引力時，固體和氣體間的介面張力會將液-固介面拉伸，換句話說，被濕潤的固體表面有較低的介面張力，因此液體會在固體表面擴張。

液體的親水性使液滴在固體表面凝結，形成較大的液滴狀，而此時液體和固體表面所形成的夾角稱為接觸角(contact angle)，接觸角為固體表面和液滴切線的夾角如圖2-6，假使接觸角小，如水滴在玻璃基板上的情形，表示液體易濕潤固體表面，但是，如果接觸角像水銀液滴在玻璃基板上那麼大，代表液體不易濕潤此表面，而濕潤張力和接觸角的關係，可以用楊格方程式(2.31)表示，氣-固與液-固介面張力的差等於氣-液介面張力乘上接觸角的餘弦函數。我們先考慮兩種極端的情形，當接觸角為0度時，表示液體能完全的濕潤於固體表面；當接觸角為180度時，代表液體完全不能濕潤於固體表面，從固─氣的介面張力觀點來看，當接觸角越小，餘弦函數cos 會越大，固-氣的介面張力也會越大，此時表示固體表面較易被濕潤，而當接觸角越大，固-氣的介面張力越低，代表越不易被濕潤。

Young's equation :

SV - LS = VL cosθ (2.31)

其中
SV surface energy between solid and vapor
LS surface energy between liquid and solid
VL surface energy between vapor and liquid
contact angle

當 = 00
代表液體能完全的濕潤於固體表面

當 = 1800
代表液體完全不能濕潤於固體表面



第三章 模擬方法

3.1分子動力學模擬流程圖

本論文所撰寫的程式主要參考幾本分子動力學相關之書籍，其中包括J. M. Haile 所著作的分子動力學模擬(Molecular Dynamics Simulation)、M. P. Allen 所著液體之電腦模擬(Computer Simulation of Liquid)、D.C.RAPAPORT所著分子動力學模擬藝術(The Art of Molecular Dynamics Simulation)以及R.J.Sadus所著液體分子模擬(Molecular Simulation of Fluids)，分子動力學模擬的電腦模擬流程圖參考表二，程式計算的流程為表三，MD模擬的流程係先建立水分子SPC/E和白金(platinum)分子的potential model，加入適當的運動方程式(Newtonian或Hamilton)，利用建立出來的model，根據第二章的理論基礎來計算分子之間的勢位能(即水分子與水分子以及水分子對白金分子)，然後利用Gear’s Predictor Correction Algorithms可以計算出分子之間運動的軌跡，利用統計熱力學進而推算出各種熱力學性質(溫度、動能、位能、徑向分佈函數…等等物理性質)，論文中所使用的繪圖軟體為PVWIN，是由日本東京大學Maruyama實驗室所提供之軟體，分子動力學模擬程式建構的基本流程可分成三個主要步驟，a.初始(initialization)、b.平衡(equilibration)、c.產生(production)，下麵詳細說明這三大步驟的流程。

3.1.1初始(initialization)

� 當選定好所需的勢位能，我們就可以開始做電腦的MD模擬，初始可以分成兩個部分：決定相關的參數和分子的初始化，決定的相關參數就是設定系統的單位以及各項參數的設定。初始化就是分子結構排列和初始條件的給定，即設定分子的初始位置和初l速度、角速度，在這一步驟中，最重要的工作就是修正系統的參數使之能符合所規範的條件，為了方便各種成分分子的探討，我們這系統的單位是利用無因次化的方式(3.3小節)，然後利用有限差分的運算方法(3.4小節)。

3.1.2平衡(equilibration)

在MD模擬中Ap算所需的性質之前，控制系統需要達到平衡狀態，若所模擬的系統為< N V T >系統，我們必須把分子的速度修正所要求的溫度；若為< N V E >系統，就把系統的總能修正到固定的能量；至於< N P T >系統，就必須把系統修正到固定的壓力和溫度，我們這裡利用平衡所跑出來的結果帶入至產值的初始值，也就是Relax的動作，使得分子的分佈情形更接近一般物理情形。本論文中所模擬的系統為< N VT >，平衡系統與非平衡系統的差異，與熱傳特性有關的，主要有二點：第一點是在< N V T >漸t峇?空間體積V、空間內分子數N、及溫度T係維持固定)，其中溫度控制(temperature control)多以調整分子動能的方式，目的在使所有的分子在固定且均勻的溫度囓韝UAF成最終漸倍聾嬪G。第二點是邊界條件的處理，在平衡系統的模擬中，多以週期性邊界(PBC)的觀念處理，目的僅在避免邊界條件的設定改變取樣的均質性，系統邊界中固體壁以鏡惜炷g(mirror reflection)處理，或以絕熱壁視之，此二方式皆無考慮固體壁對介質分子的熱傳，系統內部為平衡均勻(equilibrium and homogeneous)，即不存在有熱傳的行為。

這邊有兩種方法可以監控此系統的平衡狀態，分別是位置亂度(position disorder)和速度分佈(velocity distribution)C

(a).位置亂度(Position Disorder):

由N個分子所組成的面心結構(FCC)系統的情況，Verlet定義轉移階級參數λ，其定義為(3.1)

(3.1)

其中

(3.2)

d the length of one edge of the fcc unit cell

在開始狀態，所有的分子皆呈現FCC結構，因為所有 、 、 全為(1/2)a的倍數，所以λ=1，當系統已完全碰撞分離， 會在零上下震盪，震盪的震幅大小約為 ，因為系統的晶格碰撞分離之後，分子的位置分佈是隨機分佈，所以當系統從開始到平衡時，λ值|從1開始遞減直到0附近，然後在0附近上下震盪，而從模擬開始狀態到平衡狀態所需時間長短是和初始速度的給定有關，如果每一分子的初始速度相同，則系統需要較長的時間來達到平衡狀態，此方法不但可以監控分解之後的初始狀態，也可以監控目前的結構是否達平衡的狀態。

(b).速度分佈(Velocity Distribution)

為了監控麥斯威爾速度分佈(maxwell velocity distribution)的情況，我們須計算其速度分佈函數 ，並觀察它到平衡時的變化情峞A我們搨n一個和速度分佈有關的參數以瞭解速度分佈的情況，這裡我們用到的是Boltzmann’s H-function。

在計算H-function之前，我們必需先計算出速度分G蝻ヾA根據(3.3)，將不同速度區間(dv)的原子，利用函數的方式表示，根據卡式直角座標可以改寫成(3.4)。在MD模擬的中，我們選定一個有限的速度間格，其大小為 ，然後計算速度在 之間的原子數(3.5)，利用(3.5)可以得到某瞬態的速度分佈，求得在該瞬間的H-function (3.6)，同理可求得y方向和z方向的H-function，所以可以求得整個系統瞬態的H-function (3.7)，我們也可以利用H函數的震盪關係來監控系統的平衡狀態。

(3.3)
(3.4)
(3.5)
(3.6)
(3.7)

3.1.3產生(production)

當系統已F平衡狀態，帶入平衡步驟所跑出來的終值為產生的初始值，開始執行production，這步驟開始進行對各種性質的計算，因為平衡之後的狀態較接近自然界的物理現象，所以麂B驟再做各性質的計算C這個步驟和equilibration一樣，重複利用有限差分的方法來計算分子運動軌跡，另外這步驟過程中比平衡過程多了一個計算徑向分佈函數(radial distribution function)(3.8)的步驟，這個函數不但可以使我們瞭解分子間的分佈也可以知道分子的物理狀A(G體、氣體或是固體)，這步驟可以用來檢測系統分子的分佈是否為水分子液態的分佈情形，M為總時間步階。

(3.8)
其中
(3.9)

3.2系統架構

本論文分子動力學的模擬中，水分子是以FCC(face center cubic)結構來排定如圖3-1，水分子的各項參數可參考表一，而白金觸媒層是一層結構的FCC排列方式(圖3-2)，系統初始位置的架構如圖3-3，在周圍邊界條件的部分，我們是以週期性邊界條件處理(PBC)，我們在第二章有詳細介紹過邊界條件處理方式，前後左右為週期性邊界條件，上端為鏡面反射邊界，觸媒層部分為2704顆白金分子，白金分子彼此之間的距離為0.277nm總，而系統分子的結構在二維的結構如圖3-4，系統的長寬高皆為14.5nm的正立方體，系統溫度為350 K。

3.3物理參數無因次化

在孤立系統中MD的平衡，其中N代表系統的分子數，V代表系統的體積大小，E代表的是系統的總能，因此當我們決定N就可以經由ρ*來決定系統的體積V，為了方便各種分子的研究，系統的單位是採用無因次化的方式來設定如表四所示。經過無因次化的L-J model則形成(3.10)，無因次化的方式可以避免重複模擬的步驟，我們只需知道各種L-J分子(He、Ne、Ar、Xe…等分子)的 、 將之帶入L-J model即可求得其勢位能。根據無因次化的方法，我們可以將前面計算能量改寫成無因次化的表示式，總能量表示法為(3.11)，內能是(3.12)

(3.10)
(3.11)
(3.12)

其中
(3.13)

3.4數值方法

分子動力學的模擬之中，我們使用的模擬是假設分子為軟性球體(soft sphere)，而非硬性的球體(hard sphere)，假設分子為硬性的球體，我們所解的方程式是代數方程式(algebraic equations)，碰撞之後球體的運動為直線的軌跡，速度為定值，但是在實際的情況中，分子的物理性質叫接近軟性的球體。在MD模擬之中分子的勢位能隨著距離而變化，且是連續的變化，所以分子運動軌跡是非直線，速度也不等於常數，也就是兩分子並非瞬時的碰撞，兩分子因為排斥的關係而影響彼此一段有限的時間，所以我們需要用到數值方法來解微分方程式，這邊我們所使用的數值方法是有限差分法。

在MD的模擬之中，當我們設定好分子初始位置，我們需要設定各個分子的初始速度，我們將每一個分子隨機分配一個介於-1到+1之間的速度，因此系統內的分子速度將會介於-1到+1之間而平均分佈，接著再對系統的狀態來對每個分子的速度作修正。由於我們所模擬的系統是不受外力作用，亦即為孤立系統，所以必須將系統的總線性動量修正到零(3.14)，修正的步驟可以使系統穩定模擬。

定義系統的總能為E，下標D和A分別代表欲達到之值(Desired)和實際之值(Actual)，系統實際上總能及動能、位能之間的關係如(3.15) 若要把系統的總能修正到 ，所需要的動能為(3.16)，然後將先前隨機所給定的初始速度由方程式式(3.17)修正，將x、y、z個方向的速度修正到符合系統條件的速度值。我們使用的系統為<NVT>系統，所以要對速度作溫度上的修正，修正到欲達到的溫度 ，然後可利用動能和溫度的關係(其中水分子的動能部分包含移動所造成的動能和轉動所造成的動能)，即方程式(3.18)求得實際上的溫度 (3.19)，上式中括號 代表該值對時間的平均值，求得實際上的溫度量後，便可利用(3.20)對初始速度值作修正。

(3.14)
(3.15)
(3.16)
(3.17)
(3.18)
(3.19)
(3.20)

其中
k = 1.381x10-23 J/(molecule K)Boltzmann’s constant
vnew new velocity
vold previous velocity
I* the principle moments of inertia
W* angular velocity

3.4.1 有限差分法(Finite Difference)

泰勒定理開創了有限差分理論，使任何單變量函數都可展成羃級數，有限差分法(Finite Difference Method)的原理，是將一連續之物理區間，以有限形狀規則的網格填滿，在此離散的格點中，代入控制方程式中，並配合起始條件(initial conditions)及邊界條件(boundary conditions)，計算一格點與其附近格點間的關係。常用在分子動力學模擬的方法有，Verlet、Leap forg和Gear’s Predictor Correction，下麵一小節詳細說明Gear’s Predictor Correction。

(a).Verlet:

Verlet的方法是由前兩個時間的位置，根據(3.16)計算出下一時間的分子位置，速度則利用(3.17)求得

(3.16)
(3.17)

(b.)Leap forg：

Leap forg的方法需要知道前一時間的位置ri(t)，中間時間的速度Vi(t-0.5h)和力Fi(t)，首先由方程式(3.18)計算出Vi(t+0.5h)，然後再由(3.19)計算新時刻的位置ri(t+h)，最後由(3.20)求得速度Vi(t)

(3.18)
(3.19)
(3.20)

3.4.2 Gear’s Predictor Correction Algorithms

預測校正法是分子動力學模擬中最常用的方法，其基本思想是Taylor series展開，我們的模擬所使用的方法為Gear’s預測校正法，此方法可以求解較複雜的勢位能函數，結果也較精確，在這裡我們除了分子運動位置利用此法修正以外，quaternion和角速度也利用此方法來修正，首先我們將五階泰勒級數展開(3.21~3.26)。

(3.21)
(3.22)
(3.23)
(3.24)
(3.25)
(3.26)

修正預測的位置和其導函數，我們使用加速度 ，預測分子之間作用力的計算方式可由(3.27)求得在 的力，將此力帶入牛頓第二運動方程式(2.2)可以得到預測加速度 ，其關係是式為(3.28)，在二階的Gear’s Predictor Correction Algorithms，微分項用來修正所有預測的位置和其導函數，方程式如下(3.29~3.33) ，預測校正因數參考表五。
(3.27)
(3.28)
(3.29)
(3.30)
(3.31)
(3.32)
(3.33)
其中
(3.34)

3.4.3 Neighbor List

MD模擬中最消耗時間的部份為計算分子間作用力及勢位能，假設系統中有N個分子，電腦每次運算需要計算(1/2)N(N-1)次的分子間作用力，如果使用的分子越多對電腦資源的消耗嚴重。在第二章2.4小節，我們介紹截斷勢位能(Truncated Potential)的方法，當分子之間的距離大於截斷長度rc，則分子間的作用力及勢能是忽略不計，因此，電腦不需要去計算當分子之間距離超過rc 的分子間勢位能和相互的作用力。為了提高電腦運算的速度，我們這邊利用Neighbor List的方法來避免不必要的計算浪費，對於系統中每一個分子，我們將所有和該分子相距在特定距離rL以內的分子全部編號紀錄在同一行矩陣上，這樣我們就可以知道對該分子有影響的分子為何，而紀錄這些分子的編號的這一行矩陣我們將它稱為Neighbor List。Neighbor List的方法是利用兩個一維矩陣來儲存分子之間的位置資料，我們在碰撞一定次數之後記錄一次數據，一般而言我們設定十次，而rL的長度大小通常些微大於rC(0.3 )，Neighbor List方法需要設置兩個一行矩陣，List和 Npoint的關係如圖3-7。List用來紀錄每一個i分子以及相鄰分子j的編號，只紀錄j > i 的分子，(編號數j < i 分子的Neighbor List裡i分子已經包含在其中)，所以不需要再重覆計算， List僅紀錄每一個分子的Neighbor List，因此在資料中我們無法知道那幾些分子和i分子相鄰，所以我們利用另一行距陣紀錄那些分子是和某分子i的相鄰分子，將紀錄在List中，利用Neighbor List的方法，可以大幅的降低電腦運算所需的時間，對於模擬的速度有很大的改善，因此用在本模擬系統中如圖3-7，我們考慮系統之中有四種原子，分別是Oxygen和Hydrogen a、Hydrogen a以及Platinum四種原子，所以我們只需計算O-O、O-Ha、O-Hb、O-Pt、Ha-Ha、Ha-Hb、Ha-Pt、Hb-Hb、Hb-Pt、Pt-Pt十種分子之間的距離，根據其距離來判斷是否落在Rcut之內來決定是否計算其勢位能和力，決定之後在根據其不同原子和原子之間的勢位能來判別利用何種勢位能方程式計算其勢位能和力，利用此方法可以減少計算所需消耗的時間。

3.4.4 Quaternion for Non-Linear Molecular

由於水分子是三種原子所構成的結構，如圖2-2我們知道水是由兩個氫和一個氧所構成，因此水是屬於非線性的分子，碰撞的過程之中會有轉動(Rotation)和震動(Vibration)的物理現象，由於分子震動是在極短的時間內非常迅速的運動，因此我們忽略不計，這裡我們假設分子的鍵長是固定，不會有震動的現象，因此我們只計算移動和轉動的部分，移動的部分3.4.3小節提過，轉動的部分，根據由拉Euler運動的方程式可寫出空間中剛體角速度的方程式(3.35)、(3.36)、(3.37)，經過(3.38)可算出Euler運動方程式剛體的角度，由此式之中我們知道當sin 等於零和 時，其值發散是為奇異點，因此這裡我們使用Quaternion的方法來避開奇異點，利用四組的quaternion< Q(q0,q1,q2,q3) >參數來產生座標，其中參數滿足方程式(3.39)，我們將A重新定義為(3.41)，並滿足(3.42)，再利用Gear’s Predictor Correction的方法來修正其計算出來的quaternions參數以及角速度的一次到四次的微分值，不同於位置的修正，位置的運動方程式是屬於二階微分方程式，這裡的修正是屬於一階微分方程式，參考表六，利用代數的方法來解決奇異點的問題，並提升運算所需時間。
A=
(3.35)
(3.36)
(3.37)
(3.38)
(3.39)
(3.40)
(3.41)
(3.42)

Where

Ixx the principal moments of inertia of x axis
Iyy the principal moments of inertia of y axis
Izz the principal moments of inertia of z axis


第四章�結果與討論

本論文主要是模擬白金觸媒層和水分子反應的物理現象，根據分子動力學的基礎建立程式模擬，我們採用四種不同大小的水分子團來模擬，分別是256顆水分子、500顆水分子、864顆穭壑l和1936顆水分子，然後對於模擬得到的物理現象探討，討論接觸角的產生變化情形以及水分子擴散附著於白金觸媒層的速度和時間的關係，其結果如下：

4.1擬痕G

4.1.1 SPC/E Model
在MD的模擬中，為了確定我們所建立的水分子(SPC/E)的Model是否正確，因此我們建立好Model之寣A設定系統狀態298K，模擬液態水分子的分子動力學模擬，其中液態水的密度在298K時根據熱力學飽和液態水查表可得其密度為998 kg/m3，利用此密度，根據方程式ρ=NM /L3換算分子之間的間距，有了分子之間的間距之後，在建構水分子的結構進行模擬，N代表分子數、M代表水分子質量、L代表分子結構長度，系統模擬進行25 ps (5 ps relax和20 ps平?，每一步階為0.5fs。

我們利用兩個圖來印證水分子的Model是否合理，分別是徑向分佈函數圖以及勢位能的震盪圖，圖4-1為水分子的徑向分佈函數圖，其中其分佈的型態為水分子之分佈情形，其第一個鋒值大約在分子距離為0.25nm左右有三個氧分子的分佈情形，第二個鋒值大約在0.47nm，圖4-2為水分子的勢位能能量圖的震盪情形，其值大約在-41.2~-41.5 kJ/mol 間來回震盪，而[6]所計算的值大約為-41.4 kJ/mol，此外我們所建立的Model是以<NVT>的系統來建立，圖4-3為溫度震盪的情形，從圖中可以看到溫度的震盪趨於穩定的狀態，經過幾個簡單的驗證之後，接著說明為我們建構的模擬系統。
4.1.2 256、500、864、1936顆水壑l之系統模擬結果
計算所需的性質之前，控制系統需要達到平衡狀態，達到平衡狀態之後，我們所模擬的結果才可以更符合物理現象，當系統已達平衡狀態，這裡我們分別模擬256顆水分子、500顆水分子、864顆水分子和1936顆水分子在白金觸媒層上的反應情形，系統ensemble的選定為<NVT>，系統的邊長皆?4.5nm，一個時間步階為0.5 fs，溫度設定為350K，在給定初始狀態之後，系統先先跑200000步階(100 ps)，系統速度的修正100 ps，此步驟為equilibration，目的在使畯怐漕t統狀態更接近自然界的物理現像，100 ps之後，系統繼續模擬900 ps，這部分不做系統速度修正，此步驟為production，共1000ps（2000000步階），系統詳細的參數設定列於表六。

以下針對四種不同水分子數目大小的分子團的系統作模擬，水分子團的大小分別為256顆水分子團、500顆水分子團、864顆水分子團和1936顆水分子團，四種分子團的系統分別模擬1000 ps，圖4-4為256顆水分子團在時間0 ps的初始狀態，圖4-5是256顆水分子團在時間100 ps的模擬結果、圖4-6O256�水分子團在時間200 ps的模擬結果、圖4-7是256顆水分子團在時間300 ps的模擬結果、圖4-8是256顆水分子團在時間400 ps的模擬結果、圖4-9是256顆水分子團在伅?00 ps的模擬結果、4-10為256顆水分子團在1000 ps的模擬結果，

圖4-11為500顆水分子團在時間0 ps的初始狀態，圖4-12是500顆水分子團在時間100 ps的模擬結果、圖4-13是500顆水分子團在時間200 ps的模擬結果、圖4-14是500顆水分子團在時間300 ps的模懇痕G、圖4-15是500顆水分子團在時間400 ps的模擬結果、圖4-16是500顆水分子團在時間500 ps的模擬結果、4-17為500顆水分子團在1000 ps的模擬結果。

圖4-18為864顆水分子團在時間0 ps的初始狀態A圖4-19是864顆水分子團在時間100 ps的模擬結果、圖4-20是864顆水分子團在時間200 ps的模擬結果、圖4-21是864顆水分子團在時間300 ps的模擬結果、圖4-22是864顆水分子團在時間400 ps的模擬結果、圖4-23是864顆水分子團在時間500 ps的模擬結果、4-24為864顆水分子團在1000 ps的模擬結果，而圖4-25∼圖4-26為864顆水分子在時間1000 ps的切面圖，切面大小分別為2nm和4nm的切面圖。

圖4-27為1936顆水分子團在時間0 ps的初始狀態，圖4-28是1936顆水分子團在時間100 ps的模擬結果、圖4-30是1936顆水分子團在時間200 ps的模擬結果、圖4-31是1936顆水分子團在時間300 ps的模擬結果、圖4-32是1936顆水分子在時間400 ps的模擬結果、圖4-33是1936顆穭壑l團在時間500 ps的模擬結果、4-34為1936顆水分子團在1000 ps的模擬結果，而圖4-34∼圖4-35為1936顆水分子在時間1000 ps的切面圖，切面大小分別為2nm和4nm漱蟑措洁A而圖4-36為1936顆水分子在1000 ps的三維示意圖，圖4-37為水分子擴散在白金觸媒表面的速率圖。

4.2討論

1. 由圖4-4∼圖4-10是256顆水分子團從時間0 ps~1000 ps的模擬結果，從圖片中可以觀察到，當模擬開始隨著時〞獐W加，水分子受到白金觸媒表面的吸引力的影響，會逐漸附著於白金觸媒的表面之上，這是因為白金觸媒表面是屬於親水性(hydrophilic)，巨觀的角度來看，水滴之所以附著於白金觸媒表面之上是由於重力的關係，導致水會像白金觸媒表面的方向附著，但是在微觀的角度來看，重力已經不是影響接觸角的唯一因素，表面張力的影響已經是影響接觸角的另外一項重要因素，以奈米尺度觀點來看，重力影響幾乎微乎其微，而表面張力的影響也是其中一項因素，但是這邊我們暫時不探討這項因素，我們將討論重點集中在分子之間的作用力，因為白金觸媒表面的親水性，造成白金分子吸引水分子，形成水分子彼此之間的作用力和白金與水分子之間的吸引力相互抗衡，但是由於白金觸媒表面的吸引力大於水分子之間的庫倫力，因此水滴逐漸被附著於白金觸媒表面之上，從圖中我們可以觀察到水分子團在白金觸媒層之上隨著時間慢慢的擴散並且附著(adhesion)於白金表面之上，但是由於水分子數目不夠多，因此只有第一層(first layer)水分子薄膜的出現，觀察不到我們所想觀察的物理現象。

2. 圖4-11∼圖4-17是500顆水分子團從時間0 ps~1000 ps的模擬結果，從圖我們可以觀察到水分子團在白金觸媒層之上隨著時間慢慢的擴散並且附著於白金表面之上，在此可以觀察到附著於白金觸媒層表面的水分子數目明顯的增加，而第二層(secondary layer)水分子薄膜也已經出現，但是接觸角的資訊並不明顯，無法從此得知接觸角的任何資訊，因此我們繼續把系統的分子數目加大至864顆水分子。

3. 圖4-18∼圖4-24是864顆水分子團從時間0 ps~1000 ps的模擬結果，從圖我們可以觀察到水分子團在白金觸媒層之上隨著時間慢慢的擴散並且附著於白金表面之上，可以看到水分子在白金表面上的兩層薄膜清楚的出現，從圖4-25∼圖4-26切面洏i以看到接觸角的形成，這是因為水分子數目的增加，使得庫倫力和白金表面的吸引力大小越來越接近，造成在力平衡上的牽制而導致接觸角的形體產生，根據[13]所得到的實驗結果(表七)，水滴在白金表面上的附著情形，不同液滴大小附著所造成的角度不同，分別是22.05度、40.81度、50.71度，這裡我們取圖4-26的切面圖來量測接觸角的大小，選擇此切面的原因是因為不做切面的觀測，有可能因為角度的問題造成觀測上角度的誤差，而影響量測角度的結果，因此我們利用切面的方式來觀測模擬結果，從圖中我們利用角度規量測到 的角度為19度。

4. 我們把系統水分子數目加大至1936顆水分子，這裡我們排定的方式是屬於長條狀的水分子團，利用這種排列方式是因為排定的形狀越接近自然型態的液滴形狀，計算系統達平衡的時間也會相對減少，圖4-27∼圖4-33是1936顆水分子團從時間0 ps~1000 ps的模擬結果，從圖我們可以觀察到水分子團在白金觸媒層之上隨著時間慢慢的擴散並且附著於白金表面之上，而隨著時間的增加其接觸角的形狀也清楚的出現，隨著時間的增加接觸角的變化慢慢趨於穩定平衡狀態，變化情形會越來越小，另外圖中可以看到有一些零星的水分子散佈在空間系統中，因為系統溫度為350 k，所以判定這些水分子是屬於汽化狀態的水分子型態，伴隨著比較大的動能分佈於系統中，從圖4-34∼圖4-35切面圖可以看到接觸角的形體的大小以及形狀，接觸角的形體已經產生，因為白金觸媒的吸引力使水滴附著於表面之上，但是由於水分子數目足夠多，所以庫倫力的影響比重也比較明顯，使液滴的形狀清楚明顯，根據圖4-35切面圖，我們利用角度規量測出來接觸角 的角度為42度，將量測的角度整理於表七由於不同水分子團的大小造成量測出來的結果值不同，和實驗的結果趨勢相同，因為水分子的數目比較多，潤濕的效果明顯的增加。

5. 圖4-36為1936顆水分子的三維視點觀測圖，黃色的虛線圓圈部分是水分子團所造成的擴散範圍，藍色線代表的是擴散半徑，紫色的圓圈是小團的水分子凝聚的物理現象，圖4-37為水分子在白金觸媒層表面上的擴散速率，圖中可以觀察到水分子一開始的擴散速率比較急速，到了500 ps之後擴散速率逐漸的緩和，到達最後1000 ps時，面積幾乎不再增加，反應已達平衡狀態，圖中1936顆水分子的擴散速度(紫色曲線)開始的趨勢比其它三組快，這是因為一開始給定的初始條件比較接近自然情況，所以呈現比較和緩的擴散，實驗上[15]擴散半徑的速度大概是 ，模擬的結果大概趨近於 。


第五章�結論與建議

本論文主要探討白金觸媒和水的反應變化，透過分子動學方法所的模擬得到的結果可以讓我們瞭解微小分子之間的物理變ヾ形、同時預測腔疏云漣峖酉P變化、計算水分子擴散面積的速度，除此之外，這邊也對外來研究的方向提供幾點建議。

5.1結論

1. 由於水分子的極坋釋部A使得水分子|形成大小不一的水分子團，但是如果水分子的數目不夠多，所形成的水分子團會不夠大，而且水分子團也無法在白金觸媒層表面上形成明顯的薄膜，當水分子的數目夠多時，在模擬達平衡狀態時就會在白金觸媒層表面上形成一層薄膜(first layer)，而要形成第二層薄膜(secondary layer)，則水分子的數目至少需要500顆以上的水分子，如果要有接觸角的形成，則必須把水分子的數目加大到至少864顆水分子以上，才會有微小的形狀產生，但是並不明顯，因此要觀測到接觸角的形成則系統的水分子數目至少要有上千顆以上的數目，才能觀測到接觸角的產生以及變化情形，864顆水分子團為19度，1936顆水分子團為42度。

2. 當水分子的數目夠多，多到可以形成接觸角時，透過模擬發現接觸角的形成和系統初始排定的水分子結構形狀有密切的關係，如G我們給定的形狀越接近自然情況下水分子團的形狀，系統在模擬的時候會比較快達到穩定的狀態，如此一來可以降低模擬所需要消耗的時間。

3. 分子在白鷵眼C層上的擴散速度在500 PS前，擴散的速度比較快，但是到了500 PS之後，擴散的速度會逐漸變慢，這是因為開始時水分子團受到白金觸媒層表面的吸引力所吸引，水分子團|迅速的向外擴散於觸媒層表面，但是經過一段時間，分子之間的作用力快要達到平衡的狀態，因此擴散速度逐漸變慢，並且有第二層水分子薄膜的形成，當系統到達1000 PS時，接觸角的形體逐漸產生，系統達平衡狀態，接觸角的大小變化甚小，而水分子團擴散面積的速度趨近於 。

4. 不同的水分子團也|因為彼此之間的庫倫力l引，造成水分子團的形成，在1936顆水分子的模擬中，小團的水分子會受到大團水分子的吸引力而聚集，這是因為小團的水分子的總內能比大團水分子的總內能小而遭到吸引，而水分子附著在白金觸媒層上形成薄膜的原因是因為受到水分子和白金分子之間的吸引力而形成薄膜。

5. 接觸角的形成是因為系統已達準平衡的狀態，分子之間的內聚力和附著力達平衡，水分子與水分子之間以及水分子和白金之間的作用力以達平禳A以巨觀的角度來看，在形體上肉眼幾乎看不出有任何變化情形，但是在微觀的角度來看，分子之間的作用力不斷的在變化，水分子團持續不斷的運動(移動、轉動、震動)，接觸角的形體還是有些微的變化現H。

5.2建議

1. 架設接觸角量測的實驗，研究液滴(水、甲醇)和基材(白金、碳布)的接觸反應實驗，目前的技術在量測奈米液滴上有其困難，蒸發現象明顯量測ㄘ騿A必須等到技術成熟，再將接觸角量測反應情形和模擬結果比較，以確認模擬結果的準確性和可靠性。

2. 在燃料電池的模擬中，我們關心的問題是甲醇和碳布基材的反應情形，建議加入更複雜的甲醇分子和碳布基材來進行模擬，並且和實驗數據騆�。

3. 加入奈米液滴之相雂じ狾﹛A研究不同種類之分子凝聚(cluster)的液滴成形情形、奈米液滴蒸發及固化、相變化介面動態、表面張力等運動行為。

4. 溫度控制方法加入，phantom molecule的邊界處理方法可以達到溫度控制的目的，利用phantom molecule的方法來模擬白金觸媒層，其預測的結果與篔擘u簧的強度以及各層分子之層數(number of layers)設定值有關係，我們可視電腦效能來增加其層數，此外也可以利用constraint的方法來達到溫度控制的效果，可以[快計算所需要的時間。

5. 利用個人電腦P2.4的處理器來計算分子動力學模擬，所需消耗的時間過於冗長，因此當前如何建立平行化處理的程式，架構平行化運算的處理器為當務之急。

參 考 文 獻
1. J. M. Haile, ”Molecular Dynamics Simulation Elementary Methods”, John Wiley & Sons, 1992.

2. M. P. Allen and D. J. Tildesley, ”Computer Simulation of Liquid”, Oxford University Press, 1989.

3. D. C. RAPAPORT, ”The Art of Molecular Dynamics Simulation”, Cambridge University Press, 1995.

4. R. J. Sadus, ” Molecular Simulation of Fluids”, Elsevier Science B.V., 1999.

5. S. Maruyama, “Molecular Dynamics Method for Microscale Heat Transfer, Advances in Numerical Heat Transfer”, Taylor & Francis, New York, 2, p.p.189-226, 2002.

6. H. J. C. Berendsen, J. R. Grigera, and T. P. Straatsma, “The missing term in Effective Pair Potential ”, J. Phys. Chem., 91, p.p.6269-6271, 1987.

7. E. Spohr, “Computer simulation of the water/platinum interface ”, J. Phys. Chem., 93(16), p.p.6171-6180, 1989.

8. J. A. Nieminen, D. B. Abraham, M. Karttunen and K. Kaski, “Molecular Dynamics of a Microscopic Droplet on Solid Surface”, Phys. Rev. Lett., 69-1, p.p.124-127, 1992.

9. J. Yang, J. Koplik and J. R. Banavar, “Terraced Spreading of simple Liquid on Solid Surfaces”, Phys. Rev. A., 46-12, p.p.7738-7749

10.S. Matsumoto, S. Maruyama and H. Saruwatari, “A Molecular Dynamics Simulation of a Liquid Droplet on a Solid Surface”, ASME/JSME Thermal Engineering Conference., 2, p.p.557-562, 1995.

11.S. Maruyama, T. Kimura, T. Kurashigo and Y. Yamaguchi, “Liquid Droplet in Contact with a Solid Surface”, Micro.Thermophys Eng., 2, p.p.49-62, 1998.

12.S. Maruyama and T. Kimura, “A Molecular Dynamics Simulation of a Bubble Nucleation on Solid Surface”, Int. J. Heat & Technology, 18, p.p.69-74, 2000.

13.S. G. Kandlikar, S. Maruyama, M. E. Steinke and T. Kimura, “Measurement and Molecular Dynamics Simulation of Contact Angle of Water Droplet on a Platinum Surface”, HTD (Am. Soc. Mech. Eng.), 369, p.p.343-348, 2001.

14.S. Maruyama and T. Kimura, “Molecular dynamics Simulation of liquid Droplet on solid surface”, Proc.12th Int. Heat Transfer Conf., p.p.537-542, 2002.

15.F. Heslot, A. M. Cazabat, P. Levinson, and N. Fraysse, “Experiments on wetting on the scale of nanometers : Influence of the surface energy”, Phys. Rev. Lett., 65-5, p.p.599-602, 1990.