應用分子動力學模擬探討不同分子量對奈米尺寸液滴碰撞之影響

(1)

國立交通大學

機械工程研究所

碩士論文

Effect of Molecular Weight on Nanoscale Droplet

Collisions Using Molecular Dynamics Simulation

應用分子動力學模擬探討不同分子量對奈米尺寸液滴

碰撞之影響

研究生：王柏勝

指導教授：吳宗信

博士

(2)

應用分子動力學模擬探討不同分子量對奈米尺寸液滴碰撞之影響

Effect of Molecular Weight on Nanoscale Droplet Collisions Using

Molecular Dynamics Simulation

研究生：王柏勝

Student：Po-Sheng Wang

指導教授：吳宗信博士

Advisor：Dr. Jong-Shinn Wu

國立交通大學

機械工程學系

碩士論文

A Thesis

Submitted to Institute of Mechanical Engineering

Collage of Engineering

National Chiao Tung University

In Partial Fulfillment of the Requirements

for the degree of

Master of Science

In

Mechanical Engineering

July 2007

Hsinchu, Taiwan

中華民國九十六年七月

(3)

致謝

首先要感謝我的指導老師，吳宗信教授，不僅在每週的會議上訂正我研究的錯誤及方向，更在每天都會定時的討論我研究的進度，使我在學習的態度上有更正確的認知。同時也要感謝遠從義守大學來的江仲驊教授，中原大學的趙修武教授，成功大學的陳政宏教授願意來當我的口試委員，並對於我論文內容的批評與建議，使我的論文更佳的充實完善。同時感謝交通大學給了我求知的環境，讓我在這兩年內的碩士生涯有更多相同領域、相同的學生互相砥礪，學習。在這 700 多天的APPL實驗室裡，總是充滿了一股溫馨的氣氛。由衷地感謝李允民學長，在我對程式上有所不了解時，總會及時的給予我解答，並提醒我更多應該要注意的部份，已畢業的邵雲龍學長，增進了不少我在電腦硬體領域的相關知識。周欣芸學姐及洪捷粲學長，從你們身上我學習到對事物處理該有的守則及信條，尤其是粲哥，建立的學習方式使整個實驗室能更有規劃的朝研究之路大筆的邁進。李富利、鄭凱文、邱沅明、江明鴻及胡孟樺學長姐、還有已畢業的許國賢、陳百彥、梁偉豪及陳育進學長們，在我剛進入碩一時不時的提醒我在學校的學習課程中要注意的細節。還有我的同學們，謝昇汎、陳又寧、洪維呈及盧勁全，有了你們的協助，讓我在這兩年內學習有了更加寬廣的空間。此外，吳玟琪、鄭丞志、林士傑、蘇正勤、呂政霖、柳志良及劉育宗學弟妹，幫助我在研究時討論了許多有關未來領域方向的準備。並感謝在我準備最後論文階段時，時常關心我的朋友，科科王、走不出、祥ㄝ、國文將軍、蝗蟲及天兵顏，使得我的論文在我預想不到的情況下完成。更要提到的是博班念了七年的許祐寧，在我每次有問題時從不會告訴我正確答案，並在每次詢問時都會不斷的提醒我學習上的怠慢，但從我比較他在實驗室中及畢業後對老師的負責任的態度，更讓我對社會上形形色色的人有了更一深刻的體認。最後，我想感謝我的家人，我的父親王明亮先生，我的母親黃月珠女士，及我的弟弟王翎羽，不僅在生活上給予我全力的支柱，更在我從小至今跌跌撞撞的成績裡，不停的鼓勵我支持我，使我能一步一步的越過人生不同的里程碑。還要感謝我的女友，秝丞，這兩年有了妳時時刻刻的陪伴我協助我，使我最後完成了這篇論文，我也會以更積極的態度去面對之後的人生。王柏勝謹誌九六年七月于風城

(4)

應用分子動力學模擬探討不同分子量對奈米尺寸液滴

碰撞之影響

學生: 王柏勝指導教授: 吳宗信博士

國立交通大學機械工程學系

摘要

本文使用平行化分子動力學程式(Parallelized cellular molecular

dynamics. PCMD)來模擬兩顆在奈米尺度下(~10 nm)由氦(Helium)或

氙(Xenon)原子所構成的液滴，並採用L-J (12-6)的勢能模型來探討兩

顆液滴於真空環境時相互撞擊的行為及影響。在模擬中影響液滴碰撞

時的行為參數主要分為液滴間的相對速度、碰撞參數 (Impact

Parameter)及材質。本文模擬氦原子的相對速度範圍為250

m/s~750

m/s，氙原子的相對速度範圍為250 m/s~2250 m/s，碰撞參數則皆為

0~8.75 nm。利用可視化程式模擬觀察到的行為有：液滴結合(Direct

Coalescence)、液滴變形結合(Stretching Coalescence)、液滴拉伸破裂

(Stretching Separation)以及液滴碎裂(Shattering)。當相對速度和碰撞参

數越高時，碰撞後的液滴其破裂和旋轉的現象會較為明顯，而材質的

(5)

不同則會影響液滴碰撞時碎裂的程度，並與之前的文獻做比較來判斷

不同分子量下，奈米尺寸的液滴相互撞擊後的行為及能量變化。

(6)

Effect of Molecular Weight on Nanoscale Droplet Collisions Using

Molecular Dynamics Simulation

Student: Po-Sheng Wang Advisor: Dr. Jong-Shinn Wu

Department of Mechanical Engineering

National Chiao-Tung University

Abstract

In this thesis, Parallelized cellular molecular dynamics (PCMD) to

simulate two droplets consist of Helium or Xenon in nanoscale and

adopt the L-J (12-6) potential to discuss the behavior and effects when

two droplets collide in vacuum. In the simulation, parameters which

influence the behavior of collision primarily involve the relative

velocity between droplets, the impact parameter and the material we

use. The simulation in this context sets the relative velocity of helium

atom range from 250 m/s to 750 m/s, the relative velocity of xenon

atom range from 250 m/s to 2250 m/s, and the impact parameters all

range from 0 to 8.75 nm. By the way of visualization program “pvwin”

we can observe several behavior of simulation as follow: Direct

(7)

Shattering. The greater the relative velocity and impact parameters are,

the more obvious separation and rotation the droplets display after

collision. Furthermore differences in material will affect the degree of

shattering after collision. And we can compare the results with

literature before to study the behavior and the change of energy in

different molecular weights after the collision of droplets in the

collision of droplets in nanoscale.

(8)

List of Tables

(11)

List of Figures

Fig. 1. 1 Terminology of possible droplet-droplet collision outcome. (a)coalescence,

(b) disruption and (c) fragmentation. ...39

Fig. 1. 2 Impact parameter (b)...39

Fig. 2. 1 Cartesian frame ...40

Fig. 2. 2 Lennard-Jones (LJ) pair wise intermolecular potential ...41

Fig. 2. 3 Xenon and Helium Lennard-Jones (LJ) pair wise intermolecular potential 41 Fig. 2. 4 All pairs, (b) Cell-link, and (c) Neighbor Lists methods ...42

Fig. 2. 5 Neighbor Lists method...43

Fig. 2. 6 Cell-link + Neighbor Lists ...43

Fig. 2. 7 Periodic boundary conditions ...44

Fig. 2. 8 Proposed flow chart for parallel molecular dynamics simulation using dynamic domain decomposition. ...45

Fig. 3. 1 Head-on (b= 0) droplets pair collision initial setup on y-z plane. ...46

Fig. 3. 2 Non-head-on (ex: b= 5nm) droplets pair collision initial setup on x-y plane. 46 Fig. 3. 3 Distribution map of various regimes of Xenon droplet-collision. ...47

Fig. 3. 4 Distribution map of various regimes of Helium droplet-collision...48

Fig. 3. 5 Snapshot of Xenon droplet pair collision under vacuum, at (a) b=0, V=250m/s, (b) b=0,V=500m/s. ...49

Fig. 3. 6 Snapshot of Xenon droplet pair collision under vacuum, at (a) b=0, V=750m/s, (b) b=0, V=1250m/s. ...50

Fig. 3. 9 Snapshot of xenon droplet pair collision under vacuum, at (a) b=1.25nm, V=750m/s, (b) b=1.25nm, V=1000m/s. ...53

Fig. 3. 10 Snapshot of xenon droplet pair collision under vacuum, at (a) b=1.25nm, V=1250m/s, (b) b=1.25nm, V=1500m/s. ...54

(12)

Fig. 3. 11 Snapshot of xenon droplet pair collision under vacuum, at (a) b=1.25nm, V=1750m/s, (b) b=1.25nm, V=2000m/s. ...55 Fig. 3. 12 Snapshot of xenon droplet pair collision under vacuum, at (a) b=1.25nm,

V=2250m/s, (b) b=2.5nm, V=750m/s. ...56 Fig. 3. 13 Snapshot of xenon droplet pair collision under vacuum, at (a) b=2.5nm,

V=2250m/s, (b) b=5nm, V=250m/s. ...64 Fig. 3. 21 Snapshot of xenon droplet pair collision under vacuum, at (a) b=5nm,

V=1000m/s, (b) b=5nm, V=1250m/s.. ...66 Fig. 3. 23 Snapshot of xenon droplet pair collision under vacuum, at (a) b=5nm,

V=2000m/s, (b) b=5nm, V=2250m/s. ...68 Fig. 3. 25 Snapshot of xenon droplet pair collision under vacuum, at (a) b=6.25nm,

V=1250m/s, (b) b=6.25nm, V=1500m/s.. ...71 Fig. 3. 28 Snapshot of xenon droplet pair collision under vacuum, at (a) b=6.25nm,

(13)

Fig. 3. 30 Snapshot of xenon droplet pair collision under vacuum, at (a) b=7.5nm, V=500m/s, (b) b=7.5nm, V=750m/s. ...74 Fig. 3. 31 Snapshot of xenon droplet pair collision under vacuum, at (a) b=7.5nm,

V=1000m/s, (b) b=7.5nm, V=1250m/s.. ...75 Fig. 3. 32 Snapshot of xenon droplet pair collision under vacuum, at (a) b=7.5nm,

V=250m/s, (b) b=8.75, V=500m/s. ...78 Fig. 3. 35 Snapshot of xenon droplet pair collision under vacuum, at (a) b=8.75nm,

V=2550m/s...82 Fig. 3. 39 Snapshot of helium droplet pair collision under vacuum, at (a) b=0,

V=250m/s, (b) b=0, V=500m/s. ...83 Fig. 3. 40 Snapshot of helium droplet pair collision under vacuum, at (a) b=0,

V=750m/s...84 Fig. 3. 41 Snapshot of helium droplet pair collision under vacuum, at (a) b=1.25nm,

V=250m/s, (b) b=1.25nm, V=500m/s. ...85 Fig. 3. 42 Snapshot of helium droplet pair collision under vacuum, at (a) b=1.25nm,

V=500m/s, (b) b=3.75nm, V=750m/s. ...90 Fig. 3. 47 Snapshot of helium droplet pair collision under vacuum, at (a) b=5nm,

V=250m/s, (b) b=5nm, V=500m/s. ...91 Fig. 3. 48 Snapshot of helium droplet pair collision under vacuum, at (a) b=5nm,

(14)

Fig. 3. 49 Snapshot of helium droplet pair collision under vacuum, at (a) b=6.25nm, V=500m/s, (b) b=6.25nm, V=750m/s. ...93 Fig. 3. 50 Snapshot of helium droplet pair collision under vacuum, at (a) b=7.5nm,

V=500m/s, (b) b=8.75nm, V=750m/s. ...96 Fig. 3. 53 Snapshot of density contour and clusters size distribution of Xenon droplets collision, b=0, V=250m/s, at (a) 25ps, (b) 75ps, (c) 150ps...97 Fig. 3. 54 Snapshot of density contour and clusters size distribution of Xenon droplets collision, b=0, V=500m/s, at (a) 25ps, (b) 75ps, (c) 150ps...98 Fig. 3. 55 Snapshot of density contour and clusters size distribution of Xenon droplets collision, b=0, V=750m/s, at (a) 25ps, (b) 75ps, (c) 150ps...99 Fig. 3. 56 Snapshot of density contour and clusters size distribution of Xenon droplets collision, b=0, V=1250m/s, at (a) 25ps, (b) 75ps, (c) 150ps...100 Fig. 3. 57 Snapshot of density contour and clusters size distribution of Xenon droplets collision, b=0, V=1500m/s, at (a) 25ps, (b) 75ps, (c) 150ps...101 Fig. 3. 58 Snapshot of density contour and clusters size distribution of Xenon droplets collision, b=0, V=1750m/s, at (a) 25ps, (b) 75ps, (c) 150ps...102 Fig. 3. 59 Snapshot of density contour and clusters size distribution of Xenon droplets collision, b=0, V=2000m/s, at (a) 25ps, (b) 75ps, (c) 150ps...103 Fig. 3. 60 Snapshot of density contour and clusters size distribution of Xenon droplets collision, b=0, V=2250m/s, at (a) 25ps, (b) 75ps, (c) 150ps...104 Fig. 3. 61 Snapshot of density contour and clusters size distribution of Helium

droplets collision, b=0, V=250m/s, at (a) 25ps, (b) 75ps, (c) 150ps. ...105 Fig. 3. 62 Snapshot of density contour and clusters size distribution of Helium

droplets collision, b=0, V=500m/s, at (a) 25ps, (b) 75ps, (c) 150ps. ...106 Fig. 3. 63 Snapshot of density contour and clusters size distribution of Helium

droplets collision, b=0, V=750m/s, at (a) 25ps, (b) 75ps, (c) 150ps. ...107 Fig. 3. 64 Measurements of largest fragment of Xenon droplet pair collision, b=0,

V=250m/s, (a) Number of atoms, (b) Vibration temperature (k), (c) Rotation energy, (d) Angular momentum, respectively. ...108 Fig. 3. 65 Measurements of largest fragment of Xenon droplet pair collision, b=0,

V=500m/s, (a) Number of atoms, (b) Vibration temperature (k), (c) Rotation energy, (d) Angular momentum, respectively. ...109

(15)

Fig. 3. 66 Measurements of largest fragment of Xenon droplet pair collision, b=0, V=750m/s, (a) Number of atoms, (b) Vibration temperature (k), (c) Rotation energy, (d) Angular momentum, respectively. ... 110 Fig. 3. 67 Measurements of largest fragment of Xenon droplet pair collision, b=2.5nm, V=250m/s, (a) Number of atoms, (b) Vibration temperature (k), (c) Rotation energy, (d) Angular momentum, respectively. ... 111 Fig. 3. 68 Measurements of largest fragment of Xenon droplet pair collision, b=2.5nm, V=500m/s, (a) Number of atoms, (b) Vibration temperature (k), (c) Rotation energy, (d) Angular momentum, respectively. ... 112 Fig. 3. 69 Measurements of largest fragment of Xenon droplet pair collision, b=2.5nm, V=750m/s, (a) Number of atoms, (b) Vibration temperature (k), (c) Rotation energy, (d) Angular momentum, respectively. ... 113 Fig. 3. 70 Measurements of largest fragment of Xenon droplet pair collision, b=5nm,

V=250 m/s, (a) Number of atoms, (b) Vibration temperature (k), (c)

Rotation energy, (d) Angular momentum, respectively. ... 114 Fig. 3. 71 Measurements of largest fragment of Xenon droplet pair collision, b=5nm,

V=500m/s, (a) Number of atoms, (b) Vibration temperature (k), (c) Rotation energy, (d) Angular momentum, respectively. ... 115 Fig. 3. 72 Measurements of largest fragment of Xenon droplet pair collision, b=5nm,

V=750m/s, (a) Number of atoms, (b) Vibration temperature (k), (c) Rotation energy, (d) Angular momentum, respectively. ... 116 Fig. 3. 73 Measurements of largest fragment of Xenon droplet pair collision, b=7.5nm, V=250m/s, (a) Number of atoms, (b) Vibration temperature (k), (c) Rotation energy, (d) Angular momentum, respectively. ... 117 Fig. 3. 74 Measurements of largest fragment of Xenon droplet pair collision, b=7.5nm, V=500m/s, (a) Number of atoms, (b) Vibration temperature (k), (c) Rotation energy, (d) Angular momentum, respectively. ... 118 Fig. 3. 75 Measurements of largest fragment of Xenon droplet pair collision, b=7.5nm, V=750m/s, (a) Number of atoms, (b) Vibration temperature (k), (c) Rotation energy, (d) Angular momentum, respectively. ... 119 Fig. 3. 76 Measurements of largest fragment of Helium droplet pair collision, b=0,

V=250m/s, (a) Number of atoms, (b) Vibration temperature (k), (c) Rotation energy, (d) Angular momentum, respectively. ...120 Fig. 3. 77 Measurements of largest fragment of Helium droplet pair collision, b=0,

V=500m/s, (a) Number of atoms, (b) Vibration temperature (k), (c) Rotation energy, (d) Angular momentum, respectively. ...121

(16)

Fig. 3. 78 Measurements of largest fragment of Helium droplet pair collision, b=0, V=750m/s, (a) Number of atoms, (b) Vibration temperature (k), (c) Rotation energy, (d) Angular momentum, respectively. ...122 Fig. 3. 79 Measurements of largest fragment of Helium droplet pair collision,

b=2.5nm, V=250m/s, (a) Number of atoms, (b) Vibration temperature (k), (c) Rotation energy, (d) Angular momentum, respectively. ...123 Fig. 3. 80 Measurements of largest fragment of Helium droplet pair collision,

b=2.5nm, V=750m/s, (a) Number of atoms, (b) Vibration temperature (k), (c) Rotation energy, (d) Angular momentum, respectively. ...125 Fig. 3. 82 Measurements of largest fragment of Helium droplet pair collision, b=5nm, V=250m/s, (a) Number of atoms, (b) Vibration temperature (k), (c) Rotation energy, (d) Angular momentum, respectively. ...126 Fig. 3. 83 Measurements of largest fragment of Helium droplet pair collision, b=5nm, V=500m/s, (a) Number of atoms, (b) Vibration temperature (k), (c) Rotation energy, (d) Angular momentum, respectively. ...127 Fig. 3. 84 Measurements of largest fragment of Helium droplet pair collision, b=5nm, V=750m/s, (a) Number of atoms, (b) Vibration temperature (k), (c) Rotation energy, (d) Angular momentum, respectively. ...128 Fig. 3. 85 Measurements of largest fragment of Helium droplet pair collision,

b=2.5nm, V=250m/s, (a) Number of atoms, (b) Vibration temperature (k), (c) Rotation energy, (d) Angular momentum, respectively. ...132 Fig. 3. 89 Measurements of largest fragment of Helium droplet pair collision, b=5nm, V=250m/s, (a) Number of atoms, (b) Vibration temperature (k), (c) Rotation energy, (d) Angular momentum, respectively. ...133

(17)

Fig. 3. 90 Measurements of largest fragment of Helium droplet pair collision, b=5nm, V=500m/s, (a) Number of atoms, (b) Vibration temperature (k), (c) Rotation energy, (d) Angular momentum, respectively. ...134 Fig. 3. 91 Measurements of largest fragment of Helium droplet pair collision,

b=7.5nm, V=750m/s, (a) Number of atoms, (b) Vibration temperature (k), (c) Rotation energy, (d) Angular momentum, respectively. ...135 Fig. 3. 92 Distribution map of various regimes of Xenon, Argon, and Helium droplet

(18)

Nomenclature

D : the diameter of droplet

E : the energy

rot

E : the rotational energy

i

F : the force vector of molecular i N : number of the density

T

n

: the temperature

vib

T : the vibrational temperature

B

k : Boltzmann constant

i

m : the mass of atom i

: number of the atoms

i

r : the position vector of molecular i

i

r

JK

: the position vector of m_i

i

v

JK

: the velocity vector of m_i

ε : the strength of the interaction ρ : the density

: the length scale σ

(19)

Chapter 1 Introduction

1.1 Motivation

1.1.1 Nanoscale droplet collision dynamics

Nanoscale droplet collision dynamics plays an important role in various

technologies such as spray forming, ink-jet printing, rain drop formation, spray

scrubbing, spray cooling, surface coating, etc. In addition, the collision between two

droplets becomes the most common event in these applications. Thus, the

understanding of the fundamental collision dynamics between two droplets becomes

crucial to optimize their applications. Depending on the size of the droplets,

descriptions of the collision dynamics can be generally classified into

continuum-scale.

The normal impact of two droplets is a complicated fluid mechanics phenomenon.

The major physical processes are the conservation among energy, momentum and

angular momentum. In a collision, the droplet loses kinetic energy while it strains and

deforms. The strains lead to viscous dissipation, accounting for some conversion of

mechanical energy to heat. In advanced, surface of the droplet increases while the

original droplet breakdowns into smaller ones and surface energy increases

(20)

conversion of kinetic energy to surface energy can be viewed as a conservative

process. To continue, the droplet surface increases and surface energy increases as

well. The surface energy can be regarded as a potential energy and conservation of

kinetic energy to surface energy can be regarded as a conservative process. The

increase of surface energy during the early part of a collision will result in recoiling

and rebounding later through the conversion of surface energy back to kinetic energy.

The momentum balance occurs through a force imposed on the droplet by the wall in

a collision as the droplet loses velocity.

Restricted by the development of nanoscale technology, nowadays it is unable to

utilize the macro-vision of continuum model to explain the nanoscale phenomenon. In

order to understand the nanoscale physical mechanism, molecular dynamics (MD)

simulation is the most widespread and useful method to solve the nanoscale problem.

When discuss the behaviors of the collision dynamics between two nanoscale droplets,

the MD method can be applied to all phases of gas, liquid and solid and to interfaces

of these three phases.

1.2 Background

1.2.1 Droplet collision dynamics

(21)

Based on the understanding in the continuum scale, droplet collision can be generally

classified into several different types of collision depending on some critical

parameters, as shown in Fig. 1.1. These different types of collision process will be

introduced in the next two paragraphs.

1.2.1.1 Droplet coalescence

Droplet coalescence, which forms an integrated post-collision droplet whose mass

is equal to the sum of the mass of the pre-collision droplets, follows after the droplet

contacts. The colliding droplets coalesce when the air film thickness reaches a critical

value (~ Å, [Mackay, et al., 1963]). The droplets may coalesce temporarily or

permanently, depending on the CKE and impact parameter (b).

2

10

1.2.1.2 Disruption and fragmentation

Temporary coalescence occurs when the CKE exceeds the value for stable

coalescence and eventually results in either disruption or fragmentation. Disruption is

means that the collision product separates into the same number of droplets which

exists prior to the collision. As for fragmentation, the coalesced droplet undergoes

catastrophic break-up to form numerous small droplets.

1.2.2 Governing parameters

The character of the two droplets collision process is controlled by the Webber

(22)

meaning of impact parameter is described in the following.

1.2.2.1 Impact parameter (b)

In fig. 1.2, if we set a line (which passes trough the center of ball B) to be the x axis

and the x axis parallels the direction of incident velocity vector of ball A. The distance

from the incident path of ball A to the x axis is defined as the impact parameter (b).

1.2.2.2 Kinetic energy of collision

The kinetic energy of collision of the droplet pair with the same droplet fluid is

given by [Julius and Li, 1989]:

2 "rb" int L

p ( )

( )

2

i i i

t

E

t

E

t

m

=

+

∑

(1.2) and int

2 (3

6)

vib B

E

T

N

k

=

−

(1.3)

where R is the droplet size ratio . Where and is radius of droplet large

and droplet small, respectively. /

L S

r r r_L r_S

1.3 Literature reviews

In the past, to recognize detailed fluid mechanisms happened between the liquid

phase and the solid phase, numerous studies concerning impact of droplet on a solid

(23)

developed. Hu propose the stochastic growth of cloud droplet distributions due to

collection processes is studied using a detailed microphysical parcel model [Hu, et al.,

1998]. The results indicate that the van der Waals forces are effective in enhancing

droplet collision when the droplets are small and the distributions are narrow. Ashgriz

and Poo carried out collision experiments with water drops in the range from

micrometer size to millimeter size [Ashgriz, et al., 1990]. The two different types of

separating collisions above were identified, reflexive and stretching separating,

Harlow and Shannon are the first to simulate droplet impacting on the solid surface

[Halow, et al., 1967]. They used “marker-and-cell” (MAC) finite-difference method

to solve the fluid mass and momentum conservation equations. Tsurutani enhanced

the MAC model to include surface tension and viscosity effects, and furthermore

considered the heat transferring from a hot surface to a cold liquid droplet when it

spread on the surface [Tsuruani, et al., 1990]. Trapaga and Szekely used “volume of

fluid” (VOF) method to study impact of molten particles in a thermal spray process

[Trapaga, et al., 1990]. Fukai, Shiiba, Yamamoto, Miyatake, Poulikakos, and

Megaridis Zhao formulated a finite-element model to study the effects on the

spreading of a liquid droplet colliding with a flat surface [Fukai, et al., 1995]. Lattice

Boltzmann method excels in modeling flow problems involving multiphase materials

(24)

method for two-phase fluid flows with large density ratios and applied the method to

the simulations of binary droplet collisions for various Weber numbers and impact

parameters [Inamuro, et al., 2004]. They simulated the other types of binary droplet

collisions under certain conditions, bouncing collision for low Weber numbers and

shattering collision (Disruption or fragmentation) for high Weber numbers and

discussed the mixing processes in that two different conditions. Greenspan and Heath

studied the collision dynamics of nanometer-sized particles [Greenspan, et al., 1991].

The individual molecules were modeled as single mass particles and the

molecule–molecule interaction was described by a Lennard–Jones potential. Based on

the reviews in the above, only preliminary studies have been done in the simulation of

nanoscale droplet-droplet collision. Understanding of the droplet collision dynamics

may become important in the fast-growing nano science and technology.

1.4 Specific objectives of the proposed study

In the past, it is rare to utilize molecular dynamics method to discuss the behavior

of the liquid droplets impacting. In this work, we use molecular dynamics simulations

to compare the effect of different molecular weight on droplets, including the argon

(the molecular weight is 39.95), xenon (131.4), and helium (4.05) molecular. We use

(25)

choose Lennard-Jones (L.J.) potential to deal with interactions between xenon

droplets and helium droplets and then we apply the completed parallel MD code to

(26)

Chapter 2 Molecular Dynamics Simulation

2.1 Basic simulations of molecular dynamics

Molecular dynamic (MD) has generally used in simulating the structure of liquids,

solids, and droplets. We will introduce briefly the molecular dynamic simulation and

how to define the model applied to simulate the behavior of a droplet impacting a

solid substrate or another droplet.

First of all in a simulation, we should consider the location of all atoms. By the

location and the time derivative, all atoms can be defined their point masses in

molecular dynamic simulation. Each atom as a particle interacts with the other

particle through interaction forces derived from interaction potentials in the system,

and time evolution is governed by Newtonian mechanics. At each time step, the

acceleration of a particle used to update the position of the particle would be

calculated by Newton’s second law of motion as follows (2.1).

2 2 2 i i i i d r F mr m i d t = = (2.1)

Where r_i is the position vector of molecule as show in Fig. 2.1 . i

Since Newton’s second law is time independent or equivalently

2 2 2 i i i i d r F mr m d t = =

is invariant under time translations. Consequently, we expect there to be some

(27)

is called the Hamiltonian, H,

( N, N)

H r p =const (2.2)

Here, the momentum of molecule is defined in terms of its velocity

by . For an isolated system, total energy E is conserved, where E is equal to

the sum of kinetic energy and potential energy. Thus, for an isolated system, we

identify total energy as the Hamiltonian; then for N spherical molecules, H can be written as i p i i p = mr_i 2 1 ( , ) ( ) 2 N N N i i H r p p U r E m =

∑

+ = (2.3)

where the potential energy results from the intermolecular interactions. First

consider the total time derivative of the general Hamiltonian (2.2), ( N) U r i i i i i i dH H H H p r dt p r t ∂ ∂ = ⋅ + ⋅ + ∂ ∂

∑

∂ ∂ (2.4)

H has no explicit time dependence, then the last term on the RHS

If, as is (2.3),

of (2.4) vanishes and we are left with 0 i i i i i dH H H p r dt p r ∂ ∂ = ⋅ + ⋅ ∂ ∂

∑

_i = (2.5)

Now consider the total time derivative of the isolated-system Hamiltonian given in

(2.3), 1 0 i i i i i i dH U p p r dt m r ∂ = ⋅ + ⋅ ∂

∑

= (2.6)

On comparing (2.5) and (2.6), we find for each molecule i. i i p H r p m ∂ ₌ ₌ ∂ (2.7) and H U r r ∂ ₌ ∂ ∂ ∂ (2.8)

(28)

Substituting (2.7) into (2.5) gives 0 i i i i i i H r p r r ∂ ⋅ + ⋅ = ∂

∑

(2.9) 0 i i i i H p r r ⎛ ₊∂ ⎞_{⋅ =} ⎜ _∂ ⎟ ⎝ ⎠

∑

(2.10)

Since the velocities are all independent of one another, (2.10) can be satisfied only,

for each molecule i, we have

i i H p r ∂ _{= −} ∂ (2.11)

Equations (2.7) and (2.11) are Hamilton’s equation of motion. For a system of N

particles, (2.7) and (2.11) represent 6N first-order differential equations that are

equivalent to Newton’s 3N second-order differential equations (2.1) .

In the Newtonian view, motion is a response to an applied force. However, in the

Hamiltonian view, motion occurs in such a way as to preserve the Hamiltonian

function, where the force does not appear explicitly.

For an isolated system, the particles move in accordance with Newton’s second law,

tracing our trajectories that can be represented by time-dependent position vectors ri(t).

Similarly, we also have time-dependent momentum pi(t). At one instant, there are

positions and moment of the N particles in a 6N-dimensional hyperspace. Such a

space, called phase space, is composed of two parts: a 3N-dimensional configuration

(29)

the components of the momentum vectors pi(t). As time evolves, the points defined by

positions and momentum in the 6N phase space moves, describing a trajectory in

phase space.

2.2 Equations of motion

Potential energy represents the functions of all atoms in the system. Because of the

complexity of these functions, there is no analytical solution to the equations of

motion, therefore we must use numerical theorems in the operation. There are many

numerical algorithms has developed to calculate the integrating equations of motion.

One of them is Leapfrog method [Frenkel, et al., 1996]. Leapfrog method is the most

general method applied in MD simulation and there are some reasons for it. When we

put the leapfrog method in the integrating equations of motion, the values of

calculation will be accurate to third order. From the viewpoint of energy conservation,

it tends to be considerably better than higher-order methods. Furthermore, its

requirement of the memory storage of the computer calculation is lower than other

methods.

The introductions of Leapfrog method are as follows :

The Leap-Frog method use the velocity at 2

dt

t− and to compile the interaction

force of atoms. Use the calculated force to compute the velocity at

t

2

dt

(30)

time step, as (2.12) equation. And then apply the velocity at 2

dt

t+ to get the

location of the atom at 2

dt

t+ , as (2.13). The process of operation above is relatively

simplebecause the method we use has no need of modification.

( ) * 2 2 i i i i F t dt dt V t V t dt m ⎛ ₊ ⎞ ₌ ⎛ ₋ ⎞₊ ⎜ ⎟ ⎜ ⎟ ⎝ ⎠ ⎝ ⎠ (2.12) ( ) ( ) * ( 2 i i i dt r t+dt =r t +dt V t+ ) (2.13)

To diminish the arithmetic error, we adopt the velocity at 2

dt

t+ to solve the

location at next time step rather than the velocity at . In other words, the computed

solutions of the time average velocity between the velocity at and t would

have better accuracy than which only applied the velocity at . The calculating

processes as Fig. 2.1. t t +dt t

2.3 Potential model

It is very important to choose an appropriate potential energy model in MD

simulation when we want to approach some different realistic materials. The force

derived from the potential on an atom is due to the interaction with surroundings.

Lennard-Jones potential, one of potential methods derived from complicated

computation, will be used in the two droplets upon a solid wall model. Why we

choose Lennard-Jones potential as our main method will be explained and discussed

(31)

2.3.1 Lennard-Jones potential

Lennard-Jones potential is proposed by J.E. [Lennard Jones, 1924]. The potential

energy of a pair of atoms, and i j, locates at ri and rj as follows eq.(2.14).

( )

_ij 4 12 6 ij ij u r r r σ σ ε⎡⎢⎛ ⎞ ⎛ ⎞ ⎤⎥ = ⎜ ⎟_{⎜ ⎟} −⎜ ⎟_{⎜ ⎟} ⎢_{⎝ ⎠} _{⎝ ⎠} ⎥ ⎣ ⎦ ij rc , r ≤ (2.14)

The functional forms of the LJ and soft-sphere potentials in MD units are shown in

Fig. 2.2 (a). Where r_ij = −r_i r_j, and the parameter ε govern the strength of the

interaction and σ defines a length scale. Fig. 2.2 (b) is a relative curve of the

xenon and helium in argon Fig. 2.2 (a). The interaction repels at close range, then

attract, and is eventually cut off at some limiting separation . While the strongly

repulsive core arising from the nonbonded overlap between the electron clouds has a

rather arbitrary form, the attractive tail actually represents the van der Waals

interaction. For the interactions involve individual pairs of atoms, each pair is treated

independently, with other atoms in the neighborhood having no effect on the force

between them.

c

r

The force corresponding to u r

( )

is f = −∇u r

( )

, so the force that atom j

exerts on atom is (2.15) i 14 8 2 48 1 2 ij ij ij ij f r r r ε σ σ σ ⎡_{⎛ ⎞} _{⎛ ⎞} ⎤ ⎛ _{⎞ ⎢} =_⎜ _{⎟ ⎜ ⎟}⎜ ⎟ − ⎜ ⎟_{⎜ ⎟} ⎢ ⎥ ⎝ _{⎠ ⎝ ⎠}_⎣ _{⎝ ⎠} _⎦⎥ (2.15)

(32)

to zero, so that there is no discontinuity at r_c ; f and higher derivatives are

discontinuous, thought this has no real impact on the numerical solution.

In addiction, these units of equation are usually expressed in dimensionless format

in MD simulation. The advantage of using dimensionless format is that the equation

of motion becomes simpler because some parameters defined in the model will be

absorbed by dimensionless units. Another advantage from the viewpoint of a particle

is that transferring these units into dimensionless format can avoid some calculating

errors in computer hardware while the atoms are over-ranged.

2.4 Force computations

2.4.1 All pairs

All pairs is the simplest method to implement, but extremely inefficient when the

interaction range is relatively small compared to linear size of simulation region.

In the simulation domain, all pairs of atoms must be examined in computing process

as Fig. 2.3 (a). In fact that the amount of computation needed grows as rules

out the method for all but the smallest value of . is the number of particles.

There are two techniques which can reduce the growth rate and they will be

discussed in next two section. c r 2 ( O N ) N N 2 ( ) O N 2.4.2 Cell-link

(33)

Cell-link which can avoid most of unnecessary compute and improve the

computational effects provides a means of organizing the information about atom

positions into a form. This method is to divide the simulation region into a lattice of

small cells, and the cells should exceed in width. Then if particles are assigned to

cells on the basis of their current positions, interactions are only possible between

particles which are either in the same cells or in adjacent cells. Because of symmetry

only half the neighboring cells need to be considered. For example, a total of 14

neighboring cells must be examined in three dimensions (include the cell itself). c

r

The wraparound effects are readily incorporated into the scheme. Distinctly the

region size must be at least 4 longfor this method to be useful. There are several

ways for implementing this cell-link list method to connect the relation between

particles and cells. In the current demonstration code, it utilizes the concept of the

pointers for particles and cells. Each cell stores a particle number, which may be zero

or nonzero. Nonzero value represents a true particle number, while the zero value

represents either the last atom in the cell or an empty cell. In addition, only one array

cell List is used to represent the particles and the cells. An obvious advantage of doing

this is that we could know exactly the size of this array in advance if periodic

boundary conditions are used. Of course, there are several other methods to

implement this idea of cell-link list technique. Ideas depicted above can be clearly c

(34)

illustrated by Fig. 2.3 (b).

2.4.3 Neighbor Lists

During MD simulation, in order to consider the distance between each atom is all

in the range or not, we spend most of time in calculating the force between each

atom. Adopting Neighbor Lists can cut down the time we spend in calculating

[Lambrakos, et al., 1986 and Verlet, et al., 1967]. Neighbor List theory builds the

relationship between each atom and its surrounding atoms in a particular interval of

time steps as Fig. 2.3 (c), likes ten to twenty steps. Picking the atom as the center, is the radius of Neighbor List like Fig. 2.4 and

c

r

i

L

r r_L = + . In the particular r_c σ

interval of time steps, when we calculate the force on the atom , only use the atoms

within range to determine which atom exists in the range. And then calculate

the force of atom without calculating the distance from whole atoms between each

other at each time step in the system. For example, in the system containing

atoms, build a Neighbor List every ten time step in a simulation to establish the

relationship of distance between each two atom. The number of calculating time is

i L r r_c i N *( 1) 2 N N−

. After every atom in the region being confirmed among all atoms in

the simulation, the calculation in forces or energy of the following ten time steps only

need to estimate which atom is in the region by the atom within region. And

then, update the distributing information of itself in each atom after passing through L

r

c

(35)

the next ten time steps. When we use Neighbor Lists, we should especially notice the

choose of σ whose value should not be too less, avoiding the calculation of forces

or energy being affected by the atoms which is not within range originally but

gets into range during the ten time steps above. On the contrary, if we choose a

large value of

L

r

c

r

σ , the number of atoms in range will increase.Estimating which

atom increases its calculating frequency in the range simultaneouslyincreases the

requirementof time in calculation.

L

r

c

r

2.4.4 Cell link + Neighbor Lists

We could combine Cell link method and Neighbor List method likes Fig. 2.5 so as

to promote the performance of the simulation.

2.5 Boundary conditions

2.5.1 Periodic boundary conditions

Unless the purpose of the MD simulation is to imitate a real walls having physical

meaning, walls in simulation had better be eliminated by using Periodic boundary

conditions (PBC). The introducing of PBC is equivalent to considering an infinite

space-filling array of identical copies of simulation region. Physical meaning of

periodic boundary conditions is shown in Fig. 2.6.

(36)

Simulating in the MD system, we would like to keep the wall isothermal, so we

define a corrected layer on wall. In the study, we use the Rescaling method to modify

corrected layer. Rescaling method keeps the wall isothermal by modifying total

kinetic energy. In microcosmic point of view, the temperature is related to kinetic

energy. When we set a temperature for corrected layer, it means to set an average

kinetic energy of atoms on the corrected layer. In a word, we must keep the kinetic

energy fixed (2.16), therefore we should have a reference value. Continually, with the

use of (2.17), we compute the total kinetic energy of atoms. Finally, we start rescaling

by using (2.18) to make the total kinetic energy in the corrected layer be the same as

the reference value which we use in (2.17). 3 2 kd B d E = Nk T (2.16) 2 1 1 3 2 2 N old ka i B a i E m V Nk = =

∑

= T (2.17) * . new old kd d i i i ka a E V V V E T = = T (2.18)

2.6 Parallel molecular dynamics method

There is no doubt about that MD simulation is a useful and valuable tool. But MD

simulation is very time-consuming due to large number of time steps and possibly

(37)

the time step to be on the order of fentosecond. Many hundreds of thousand or even

millions of time steps are needed to simulate a nanosecond in “real” time scale. In

addition, up to hundreds of thousand or millions of atoms are needed in the MD

simulation, even for a system size in the nanometer scale. In the past, there have been

considerable effort [Plimpton, 1995] that concentrated on parallelizing MD simulation

on the memory-distributed machine by taking the inherently parallelism e.g.,

[Boghosian, et al., 1990] and [Fox, et al., 1998]. Generally, parallel implementation

of the MD method can be divided into three categories, including the atom

decomposition, the force decomposition and the spatial decomposition among

processors.

2.6.1 Atomic-decomposition algorithm

In the atom decomposition method, each processor, which owns nearly the same

number of atoms as other processors and in which atoms are not necessarily

geometrically nearby, integrates the Newton’s equation for all atoms and moves the

atoms of their owns. However, this method requires global communication at each

time step, which becomes unacceptably expensive as compared with the “useful” MD

computation when the number of atoms increases to a certain amount, since each

processor has to know all information (position and velocities) of all atoms at each

(38)

atoms in the system that is independent of the number of the processors, P. Thus, the

atom decomposition method is generally suitable for small-scale problem

2.6.2 Force-decomposition algorithm

In the force decomposition method, it is based on a block-decomposition of the

force matrix rather than a row-wise decomposition in the atom-decomposition method.

It improves the O(N) scaling g to be O(N / P) . It generally performs much better that

the atom decomposition method; however, there exists some disadvantages. First, the

number of processors has to be square of an integer. Second, load imbalance may

become an issue. From previous experience [Plimpton, 1995], it is suitable for small-

and intermediate-size problems.

2.6.3 Spatial-decomposition algorithm

In the spatially static domain decomposition method, simulation domains are

physically divided and distributed among processors. This method so far represents

the best parallel algorithm for large-scale problem in MD simulation for short-ranged

interaction [Karypi et. al., 1998]; however, it only works well for a system, in which

the atoms move only a very short distance during simulation or possibly distribute

uniformly in space. MD simulation of solids represents one of the typical examples.

In contrast, if the distribution of the atoms tends to vary very often in the

(39)

during simulation, which detriments the parallel performance. Thus, a parallel MD

method capable of adaptive domain decomposition may represent a better solution for

resolving this difficulty.

2.6.4 PCMD (Parallel Cellular Molecular Dynamics) algorithm

A new parallel algorithm for MD simulation, named parallel cellular molecular

dynamics (PCMD), is developed by MuST (Multiscale Science & Technology)

laboratory in NCTU in Taiwan, employing dynamic domain decomposition to address

the issue of load imbalance among processors in the spatially static

domain-decomposition method. We focus on developing a parallel MD method using

dynamic domain decomposition by taking advantage of the existing link-cells as

mentioned earlier. In this proposed method, not only are the cells used to reduce the

cost for building up neighbor list, but also are used to serve as the basic partitioning

units. Similar idea has been applied in the parallel implementation of direct simulation

Monte Carlo (DSMC) method [Nicol et. al., 1988], which is a particle simulation

technique often used in rarefied gas dynamics. Note that in the following IPB stands

for interprocessor boundary. General procedures asFig. 2.7 in sequence include:

1. initialize the positions and velocities of all atoms and equally distribute the atoms

among processors;

(40)

followed by communicate cell/atom data between processors, renumber the local

cell and atom numbers, and update the neighbor list for each atom due to the data

migration;

3. Receive positions and velocities of other atoms in the neighbor list for all cells

near the IPB;

4. Compute force for all atoms;

5. Send force data to other atoms in the neighbor list for all cells near the IPB;

6. Integrate the acceleration to update positions and velocities for all atoms;

7. Apply boundary conditions to correct the particle positions if necessary;

8. Check if preset total runtime is exceeded. If exceeded, then output the data and

stop the simulation. If not, check if it is necessary to rebuild the neighbor list of all

atoms using the most update atom information.

9. If it is necessary to rebuild the neighbor list (N=8 in the current study), then

communicate atom data near the IPB and repeated the steps 2-8. If not necessary,

then repeat steps 3-8.

In the above, in addition to the necessary data communication when atoms cross the

IPB and particle/cell data near the IPB, there are two more important steps in the

proposed parallel MD method as compared with the serial MD implementation. One

(41)

(42)

Chapter 3 Simulation of xnone and helium droplet-droplet

collision dynamics

3.1 Simulation conditions

This content mainly concerned the simulation of the behavior of xenon or helium

droplet pair collision. In order to compare the influence of crucial parameter, impact

parameter and the relative velocity between two droplets in the collision, the diameter

of single droplet was set to be similar about 10.5nm in both xenon and helium. In

detail, in the composition of xenon droplet, we arranged the L-J potential atom in

FCC crystal structure at first, which density was 0.2 and maintained equilibrium at

166.73K, and then used the PCMD to simulate the variation in a hundred thousand

time steps, and eventually obtained a xenon droplet which contained 2800 atoms. By

the same approach, the density of helium crystal structure was 1.0 and it maintained

equilibrium at 0.55K. After using PCMD in a hundred thousand time steps, we would

obtain a droplet containing 2700 atoms. To establish the simulation model of droplet

pair collision, we copied the relative coordinates of droplets and adjusted the relative

velocities between them. In addition, if the test case was head-on collision, the system

dimension should be 200σ *500σ *500σ (Fig. 3.1). For detailed observation of the

(43)

600σ *300σ *300σ (Fig. 3.2). 3.1.1 Test conditions

When simulating a collision of xenon droplet, the impact parameter was 0 and the

relative velocity ranged from 1250m/s to 2250m/s in a head-on condition, but in a

non-head-on condition, the impact parameter was 1.25nm to 8.75nm and the relative

velocity ranged from 250 m/s to 2250m/s. In simulating a collision of helium droplets,

the impact parameter was 0, but in a non-head-on condition, the impact parameter was

1.25nm to 8.75nm. And both of two conditions had a relative velocity ranging from

250m/s to 750m/s. In all of thee cases above, the original distance between droplets

were about 140Å which was shown in Fig. 3.3 to Fig. 3.4.

To identify the droplet pair behavior easily and in objectivity, we classify each

collision behavior in following manner [Svanberg, et al., 1998]:

i. Coalescence: The largest fragment contains > 80% of the molecules.

ii. Stretching Coalescence: The largest fragment contains > 80% of the molecules,

the droplet shape transform form a ball to a rotational bar, and never breakup.

iii. Stretching Separation: The largest fragment contains < 60% of the molecules,

while the sum of the two largest fragments consists of > 90% of the molecules.

iv. Shattering: The largest fragment contains < 40%of the molecules, and the sum of

(44)

v. Bounce: The droplet pair never touch each other in simulation, while the two

largest fragments consists of > 80% of the molecules.

3.2 Results and discussion

3.2.1 The Xenon droplets collision

By the way of visualization program “pvwin”, we could observe the behavior of

collision in 500ps in head-on cases (Fig. 3.5 to Fig. 3.8) and nonhead-on cases (Fig.

3.9 to Fig. 3.38). In head-on cases, the relative velocity was set from 1250m/s to

2250m/s and we found that the impacting time became earlier and the ”disk” resulting

from collision became larger while the relative velocity increased. On account of the

high molecular weight of xenon, the attraction of L-J potential was strong which

caused the ”disk” coalesce to a larger droplet. When the relative velocity reached

2250m/s, the ”disk” appeared shattering resulting in lots of small droplets.

However, in the nonhead-on cases, the relative velocity was set from 250m/s to

2250m/s. Through Fig. 3.9 to Fig. 3.12 (a), with the increasing of the relative velocity,

we observed several different events: When the impact parameter was 1.25nm and the

relative velocity was 1500m/s, it resulted in stretching coalescence, but when the

relative velocity reached 2000m/s, although shattering happened after collision, it

(45)

to Fig. 3.20 (a), the same as above, direct coalescence, stretching coalescence,

stretching separation, and shattering occurred in succession with the increase of

relative velocity. Besides, the greater the impact parameter was, the earlier the

stretching coalescence happened. Nevertheless, if the impact parameter was greater

than 5nm (Fig. 3.20 (b) to Fig. 3.38), shattering would not happen in 2250m/s relative

velocity. When the impact parameter increased to 8.75nm, coalescence still occurred

in 250m/s relative velocity and stretching coalescence still occurred in 500m/s relative

velocity.

3.2.2 The Helium droplets collision

Due to the tendency not to coalesce of the low-molecular weight atom, we used up

to 27000 helium droplets in low temperature MD simulation. In head-on cases (Fig.

3.39 to Fig. 3.40), when the relative velocity was 250m/s, droplet area increased and

then coalesce to a larger droplet as time went by. If the relative velocity was 500m/s

and 750m/s, a phenomenon like ”net” occurred and then shattered in to pieces consist

of many small droplets as time went by. Besides, if the impact velocity was higher,

the droplets which composed small shattering became smaller and dispersed evenly.

In nonhead-on cases (Fig. 3.41 to Fig. 3.52), when the impact parameter was 1.25nm

and the relative velocity was 250m/s, the droplet occurred coalescence and rotated

(46)

became stick-like at first and then shattered in to pieces consist of many small

droplets as time went by. This kind of phenomenon called stretching coalescence

would stretch and rotate when the impact parameter increased to 3.75nm and 5nm and

the relative velocity was 250m/s. In advanced, if the relative velocity reached to

500m/s and 750m/s, the stick-like phenomenon would become horizontal. As the

impact parameter was 6.25nm and the relative velocity was 250m/s and 500m/s,

stretching separation occurred, but at 750m/s shattering remained occurred. Finally,

when the impact parameter was 7.5nm and 8.75nm, the relative velocity was 250m/s,

500m/s and 750m/s, stretching separation occurred.

3.2.3 Data analysis

Fig. 3.53 to Fig. 3.60 was a display which snapshot the density contour and clusters

size distribution at 25ps, 50ps and 150ps from y-z plane of xenon droplet pair

collision in head-on cases. At relative velocity 250m/s, the droplet distribution had no

obvious changes after collision and only few pieces formed. Until the relative velocity

reached 750m/s, while the droplets collided and the coalesced, the surface area tended

to be larger and density of the area close to the droplet surface was smaller, but no

more pieces formed nevertheless. When the relative velocity was 1250m/s to 1750m/s,

the number of the pieces tended to increase as the velocity increases and we could

(47)

droplets increased and the density of the center of the droplets decreased after

collision. But at 150ps, the droplet which stretched because of collision shrank to a

smaller one and its density of the center increased also. When the relative velocity

was 2000m/s, although a phenomenon like net occurred, the number of pieces

increased substantially. But at 150ps, we could observe that the droplets coalesced

gradually and some of the pieces became larger ones because of coalescence. When

the relative velocity was 2250m/s, the pieces formed by the droplets after collision

distributed evenly along the center of collision and the size of fragment was more

average at 150ps than at 75ps. Fig. 3.61 to Fig. 3.63 was a display which snapshot the

density contour and clusters size distribution at 25ps, 50ps and 150ps from y-z plane

of helium droplet pair collision in head-on cases. At relative velocity 250m/s, although

coalescence occurred, the helium droplet had more average density from the center to

the surface in density contour at 150ps compared with the xenon droplet coalescence

and also had larger entirely surface area, but fewer fragments formed. However, when

the relative velocity was 500m/s and 750m/s, the surface area of droplet was

apparently larger as relative velocity was higher at 75ps and fragments formed more

easily. In Density contour at 150ps, we could observe that along the center of the

droplet, the size of the fragments became smaller and the number of the fragments

(48)

Fig. 3.64 to Fig. 3.75 were the data analysis of xenon droplet pair collision at

relative velocity 250m/s, 500m/s , and 750m/s and impact parameter 0, 25nm, 50nm,

and 75nm[Liu, et al., 1997 and marcus, et al., 1998]. Fig. 3.64 (a) was the atom variation of the biggest fragment from droplet and through this figure we could find

that the atom number was raise by the droplet coalescence. Fig. 3.64 (b) was the

variation of temperature (K) and through it we could find that at 25ps, it was the

collision between xenon atoms that made the temperature increase, but in 25ps to

100ps, it was coalescence that made the temperature descend substantially, however

the evaporation made the temperature tend to be raise. Fig. 3.64 (c) was the rotation

energy change of the biggest droplet as the time increased. And Fig. 3.64 (d) was the

angular momentum distribution of different direction as the time increased. Compared

Fig. 3.64 to Fig. 3.66, we could find out that as the velocity increased, the temperature

was raised apparently at collision and both the rotation energy and angular momentum

of every direction increased. Fig. 3.67 to Fig. 3.69 was the data analysis at impact

parameter 25nm. Both of them were coalescence cases, therefore the atom number of

the biggest droplet tended to decrease because of evaporation and the temperature

increased as the relative velocity of collision increased. Moreover, on account of

non-head cases, the rotation energy increased substantially after collision and the

(49)

Compared Fig. 3.70 to Fig. 3.75, we could find that the rotation energy and the

angular momentum of y-direction also rose as the impact parameter increased. And

Fig. 3.76 shown that at 275ps, the rotation energy descended on account of stretching

separation and followed by the substantially decreasing of the angular momentum of

y-direction. Fig 3.76 to Fig 3.87 shown the data analysis of helium droplet pair

collision when the relative velocities were 250m/s, 500m/s, and 750m/s and the impact

parameters were 0, 25nm,50nm, and 75nm. In head-on case, the variation tendency of

temperature, rotation energy and angular momentum result in coalescence were the

same as the variation tendency of xenon droplet coalescence at relative velocity

250m/s. However, at 500m/s and 750m/s, shattering occurred result in the decrease of

the atom number of the biggest droplet and the decrease of temperature, and that

rotation energy and angular momentum had diversified variation result from the

collision between the biggest droplet and other fragment. In non-head on cases, Fig.

3.79 and Fig. 3.82 had the same variation tendency as stretching coalescence cases of

xenon droplet. Fig. 3.84 to Fig. 3.87 shown the stretching separation of droplets

including the quickly decrease of temperature after the droplets separated, the

decrease of rotation energy and the trend which the angular momentum became gentle

on y-direction. On the contrary, in the other cases which result coalescence, when the

應用分子動力學模擬探討不同分子量對奈米尺寸液滴碰撞之影響

國立交通大學

機械工程研究所

碩士論文

Effect of Molecular Weight on Nanoscale Droplet

Collisions Using Molecular Dynamics Simulation

應用分子動力學模擬探討不同分子量對奈米尺寸液滴

碰撞之影響

研究生：王柏勝

指導教授：吳宗信

博士

應用分子動力學模擬探討不同分子量對奈米尺寸液滴碰撞之影響

Effect of Molecular Weight on Nanoscale Droplet Collisions Using

Molecular Dynamics Simulation

研 究 生：王柏勝

Student：Po-Sheng Wang

指導教授：吳宗信 博士

Advisor：Dr. Jong-Shinn Wu

國 立 交 通 大 學

機 械 工 程 學 系

碩 士 論 文

A Thesis

Submitted to Institute of Mechanical Engineering

Collage of Engineering

National Chiao Tung University

In Partial Fulfillment of the Requirements

for the degree of

Master of Science

In

Mechanical Engineering

July 2007

Hsinchu, Taiwan

中華民國九十六年七月

致謝

應用分子動力學模擬探討不同分子量對奈米尺寸液滴

碰撞之影響

學生: 王柏勝 指導教授: 吳宗信 博士

國立交通大學機械工程學系

摘要

本文使用平行化分子動力學程式(Parallelized cellular molecular

dynamics. PCMD)來模擬兩顆在奈米尺度下(~10 nm)由氦(Helium)或

氙(Xenon)原子所構成的液滴，並採用L-J (12-6)的勢能模型來探討兩

顆液滴於真空環境時相互撞擊的行為及影響。在模擬中影響液滴碰撞

時 的 行 為 參 數 主 要 分 為 液 滴 間 的 相 對 速 度 、 碰 撞 參 數 (Impact

Parameter)及材質。本文模擬氦原子的相對速度範圍為250

m/s~750

m/s，氙原子的相對速度範圍為250 m/s~2250 m/s，碰撞參數則皆為

0~8.75 nm。利用可視化程式模擬觀察到的行為有：液滴結合(Direct

Coalescence)、液滴變形結合(Stretching Coalescence)、液滴拉伸破裂

(Stretching Separation)以及液滴碎裂(Shattering)。當相對速度和碰撞参

數越高時，碰撞後的液滴其破裂和旋轉的現象會較為明顯，而材質的

不同則會影響液滴碰撞時碎裂的程度，並與之前的文獻做比較來判斷

不同分子量下，奈米尺寸的液滴相互撞擊後的行為及能量變化。

Effect of Molecular Weight on Nanoscale Droplet Collisions Using

Molecular Dynamics Simulation

Student: Po-Sheng Wang Advisor: Dr. Jong-Shinn Wu

Department of Mechanical Engineering

National Chiao-Tung University

Abstract

In this thesis, Parallelized cellular molecular dynamics (PCMD) to

simulate two droplets consist of Helium or Xenon in nanoscale and

adopt the L-J (12-6) potential to discuss the behavior and effects when

two droplets collide in vacuum. In the simulation, parameters which

influence the behavior of collision primarily involve the relative

velocity between droplets, the impact parameter and the material we

use. The simulation in this context sets the relative velocity of helium

atom range from 250 m/s to 750 m/s, the relative velocity of xenon

atom range from 250 m/s to 2250 m/s, and the impact parameters all

range from 0 to 8.75 nm. By the way of visualization program “pvwin”

we can observe several behavior of simulation as follow: Direct

Shattering. The greater the relative velocity and impact parameters are,

the more obvious separation and rotation the droplets display after

collision. Furthermore differences in material will affect the degree of

shattering after collision. And we can compare the results with

literature before to study the behavior and the change of energy in

different molecular weights after the collision of droplets in the

collision of droplets in nanoscale.

Table of Contents

List of Tables

List of Figures

研究生：王柏勝

指導教授：吳宗信博士

國立交通大學

機械工程學系

碩士論文

學生: 王柏勝指導教授: 吳宗信博士

時的行為參數主要分為液滴間的相對速度、碰撞參數 (Impact