Close
About
FAQ
Home
Collections
Login
USC Login
Register
0
Selected
Invert selection
Deselect all
Deselect all
Click here to refresh results
Click here to refresh results
USC
/
Digital Library
/
University of Southern California Dissertations and Theses
/
Shock-induced nanobubble collapse and its applications
(USC Thesis Other)
Shock-induced nanobubble collapse and its applications
PDF
Download
Share
Open document
Flip pages
Contact Us
Contact Us
Copy asset link
Request this asset
Transcript (if available)
Content
SHOCK-INDUCED NANOBUBBLE COLLAPSE AND ITS APPLICATIONS
by
Mohammad Hossein Vedadi
A Dissertation Presented to the
FACULTY OF THE USC GRADUATE SCHOOL
UNIVERSITY OF SOUTHERN CALIFORNIA
In Partial Fulfillment of the
Requirements for the Degree
DOCTOR OF PHILOSOPHY
(PHYSICS)
May 2013
Copyright 2013 Mohammad Hossein Vedadi
ii
Acknowledgements
I am truly indebted to my adviser Professor Stephan Haas for his support from the
time of application until now. Also, I would like to thank my committee members
Professor Werner Däppen, Professor Moh El-Naggar, Professor Vitaly Kresin, and
Professor Jin Ma. Professor Däppen was always there to help in times of need. In
Addition, I would like to express my deepest appreciation to Professor Aiichiro Nakano
and Dr. Ken-ichi Nomura for their kind help with my research.
iii
Table of Contents
List of Tables ..................................................................................................................... iv
List of Figures..................................................................................................................... v
Abstract.............................................................................................................................. ix
1 Introduction................................................................................................................. 1
1.1 Outlook ........................................................................................................... 1
1.2 Shock-induced Collapse of Nanobubbles ....................................................... 2
1.3 Applications of Shock-induced Collapse of Nanobubbles ............................. 4
2 Reactive Force Field Molecular Dynamics Simulations on Parallel Computers........ 8
2.1 Molecular Dynamics Simulations................................................................... 8
2.1.1 Time Integration.............................................................................................. 8
2.1.2 Periodic Boundary Conditions and Minimum Image Convention ............... 11
2.1.3 Potential Cutoff Distance and Linked-list Cell Method ............................... 14
2.1.4 Shifted and Shifted-Force Potential for Pairwise Interactions...................... 16
2.1.5 Molecular Dynamics in Various Ensembles................................................. 17
2.1.6 Physical Quantities........................................................................................ 21
2.2 ReaxFF: A Reactive Force Field................................................................... 24
2.2.1 Bond Order Calculation ................................................................................ 24
2.2.2 Charge Calculation........................................................................................ 26
2.2.3 Potential Energy............................................................................................ 28
2.3 Parallel Molecular Dynamics Simulations.................................................... 39
2.3.1 Spatial Decomposition.................................................................................. 41
2.3.2 Communication Pattern ................................................................................ 42
2.3.3 Performance of Parallel Programs ................................................................ 45
3 Structure and Dynamics of Shock-induced Nanobubble Collapse in Water ............ 47
3.1 Simulation Setup........................................................................................... 47
3.2 Structure of Water under Shock.................................................................... 50
3.3 Dynamics of Shock-induced Nanobubble Collapse ..................................... 53
4 Poration of Lipid Bilayers by Shock-induced Nanobubble Collapse ....................... 62
4.1 Simulation Setup........................................................................................... 62
4.2 Effects of Nanobubble Collapse on Lipid Bilayers ...................................... 65
4.3 Poration and Recovery.................................................................................. 71
5 Mechano-chemical Pathways to H
2
O and CO
2
Splitting.......................................... 76
5.1 Simulation Setup........................................................................................... 76
5.2 Ionization of Water ....................................................................................... 78
5.3 Mechano-chemical Reactions between H
2
O and CO
2
.................................. 82
6 Conclusions............................................................................................................... 89
References......................................................................................................................... 91
iv
List of Tables
Table 2.1: Velocity Verlet Algorithm............................................................................... 11
Table 2.2: Pseudocode for constructing head and lscl arrays........................................... 16
Table 2.3: Parameters in bond order functions ................................................................. 26
Table 2.4: Classification of interactions in ReaxFF.......................................................... 29
Table 2.5: Parameters in 1-body energy functions ........................................................... 31
Table 2.6: Parameters in 2-body energy functions ........................................................... 34
Table 2.7: Parameters in 3-body energy functions ........................................................... 37
Table 2.8: Parameters in 4-body energy functions ........................................................... 39
v
List of Figures
Fig. 1.1: Three stages of shock-induced collapse of a microbubble are shown
schematically in the left panel. When the shock front strikes on the bubble, it shrinks
because of outside overpressure, and center of the bubble accelerates upward. At this
time, liquid surrounding the bubble enter into it and forms a jet, which is focused
towards the distal point of the bubble. Finally, the jet breaks and leaves behind a tiny
bubble. Photographs of the bubble before and after passage of the shock wave are
shown in the right panel. Picture courtesy of (Ohl and Ikink 2003)........................... 3
Fig. 2.1: Schematic illustration of periodic boundary condition in two dimensions. The
MD cell is marked by gray, and its nearest-neighbor images are shown. If a particle
exits the MD box (blue particle), its image particle (red particle) enters the MD box
................................................................................................................................... 13
Fig. 2.2: Schematic illustrationof the minimum image convention in two dimensions.
Particle i does not interact with particles shown in green, and only interacts with the
particle shown in blue. .............................................................................................. 14
Fig. 2.3: Data structure for linked list method. Here, “E” stands for Empty. Picture
courtesy of Aiichiro Nakano..................................................................................... 16
Fig. 2.4: Schematic illustration of atomic configurations in energy terms: (a) 1-body, (b)
2-body, (c) 3-body, (d) 4-body, (e) hydrogen-bonding and (f) non-covalent
iterations, respectively. A red sphere represents the position of the primary atom in
each energy term. Open bars represent covalent bonds, while dotted lines are non-
covalent bonds. Picture courtesy of (Nomura, Kalia et al. 2008)............................. 29
Fig. 2.5: Schematic illustration of correspondence between vid and sid. L
x
, L
y
, and L
z
denote lengths of the MD box in x, y, and z directions, respectively. These lengths
will be used to determine the vid of the processor to which an atom belongs and then
the relation between serial and vector IDs will determine the sid of the processor.
Picture courtesy of Aiichiro Nakano......................................................................... 42
Fig. 2.6: Schematic of communication pattern in two dimensions. First, each domain is
extended from left and right. Then, the extended domain is extended again from up
and down directions. Picture courtesy of Aiichiro Nakano. ..................................... 44
Fig. 3.1: Schematic of the simulation cell. The grey plate is the momentum mirror. ...... 49
Fig. 3.2: Experimental (blue crosses) and MD (red circles) Hugoniot compression curves
for particle velocities between 0.4 and 3 km/s. The inset shows the simulation cell
and the momentum mirror (grey plate)..................................................................... 51
vi
Fig. 3.3: Oxygen-oxygen pair-correlation functions, g(r), for ambient and shocked water
with and without a nanobubble as a function of distance at u
p
= 3 km/s. The
nanobubble diameter is 10 nm. ................................................................................. 52
Fig. 3.4: Snapshot of a water molecule cluster that forms an ice-VII-like structure in the
compressed region at u
p
= 1 km/s. Red and white spheres represent oxygen and
hydrogen atoms, respectively. Dotted lines indicate the directions of nearest
neighbor oxygen atoms (within 3Å from the central oxygen atom)......................... 53
Fig. 3.5: (a) Normalized nanobubble volume, Ω(t), versus time, t, for nanobubbles of
initial sizes D = 6, 8 and 10 nm. Insets show snapshots of the largest nanobubble
before the shock wave strikes (t = 1.7 ps) and before the nanobubble collapses (t =
3.3 ps). (b) and (c) show velocity vector fields in a cross section of the system just
before (t = 3.3 ps) and soon after (t = 3.9 ps) the largest nanobubble collapses. The
velocities are measured relative to the reference frame of momentum mirror. In (b)
the region inside the white dash lines is enlarged to show a focused nanojet in the
inset; in (c) the lateral flow indicates the onset of nanojet disintegration. ............... 56
Fig. 3.6: Time evolution of the shape of the nanobubble of diameter D = 10 nm at u
p
= 3
km/s........................................................................................................................... 57
Fig. 3.7: Rotational (a) and vibrational (b) temperature profile along the direction of
shock propagation—x. Rise in vibrational temperature lacks the one in rotational
temperature. The rise in rotational temperature occurs at t = 2.9 ps while the rise in
vibrational temperature takes place at t = 3.4 ps....................................................... 58
Fig. 3.8: Effects of particle velocity and nanobubble size on nanobubble collapse.
Snapshots in panels (a)-(c) are taken for u
p
= 3 km/s and nanobubble diameter D = 6,
8 and 10 nm, respectively. Snapshots in panels (d)-(f) are for a fixed nanobubble
size (10 nm) and u
p
= 2.5, 3.0 and 3.5 km/s, respectively. The size of the high-speed
region (ellipsoidal red regions) increases more with an increase in D than with u
p
.
All snapshots are taken immediately after the nanobubbles collapse....................... 59
Fig. 3.9: [(a)-(c)] show the primary shock wave at t = 1.6 ps, 2.9 ps, and 3.3 ps,
respectively. Here, the particle velocity, u
p
= 3 km/s and D = 10 nm. [(d)-(e)] show
the water hammer shock at t = 3.5 ps, 3.7 ps, and 3.9 ps, respectively. The color
code corresponds to pressure and shock waves are identified by pressure
discontinuities. The length scale bar for these panels is shown in (a). ..................... 61
Fig. 4.1: Shock velocity versus particle velocity. The simulation results for SPC water are
in good agreement with experimental data. The inset shows the setup for shock
simulation. The gray plate is the momentum mirror................................................. 64
vii
Fig. 4.2: Schematic of the MD box containing a nanobubble in water and a lipid bilayer.
The gray plate is the momentum mirror. .................................................................. 65
Fig. 4.3: Snapshots of velocity profile for the system with D = 40 nm, T
i
= 300 K and u
p
= 0.7 km/s. Arrows show the direction of average molecular velocities and the
velocity magnitudes are color-coded. Panel (a) shows a nanojet in the system at t =
20 ps. The white vertical region is the bilayer. Panel (b) shows a spreading flow at t
= 24 ps resulting from the impact of the nanojet on the lipid bilayer....................... 67
Fig. 4.4: Average lateral velocity of water molecules in the vicinity (1nm) of a lipid
bilayer versus the distance from the center of the bilayer for a nanobubble of
diameter 40 nm. The inset shows the results for a nanobubble of diameter 10 nm. In
both systems u
p
= 0.7 km/s and T
i
=300 K................................................................ 68
Fig. 4.5: (a) and (b) are snapshots of the density of water at t = 20 and 28 ps. Here D = 40
nm, T
i
= 300 K and u
p
= 0.7 km/s. The central blue region is the lipid bilayer. Panel
(a) shows the nanojet traveling towards the distal side of the nanobubble. Panel (b)
shows the deformed bilayer and water hammer shock. ............................................ 69
Fig. 4.6: Temporal changes in the normalized order parameter for the lipid bilayer region
directly impacted by the nanojet. At t = 16 ps, the density of water molecules around
the bilayer starts to increase, and hence the bilayer begins to show increasing
disorder. .................................................................................................................... 71
Fig. 4.7: Poration of lipid bilayers by collapsed nanobubbles. Here u
p
= 0.7 km/s and D =
40 nm. The bilayer was initially in the gel phase at T
i
= 300 K (a) and in the liquid
phase at T
i
= 323 K (b).............................................................................................. 72
Fig. 4.8: Passage of a water molecule (white dot) across the lipid bilayer for the system
with D = 10 nm, T
i
= 323 K, and u
p
= 1.0 km/s. The water molecule, initially in the
red hydrophilic region on the left side of the bilayer, crosses the hydrophobic region
(light green) and ends up on the right side (red) of the bilayer................................. 74
Fig. 4.9: Deformed (a) and relaxed (b) lipid bilayer configurations. After relaxation, the
deformation in the bilayer disappears....................................................................... 75
Fig. 5.1: Schematic of the primary system. All atoms are given a velocity u
p
in the
negative x direction towards a momentum mirror, and the x component of their
velocities is reversed when they pass the momentum mirror. The nanobubble with
diameter 10 nm contains 532 CO
2
molecules........................................................... 77
Fig. 5.2: The number of ions as a function of time in the reference system (a) and in the
primary system (b). ................................................................................................... 79
viii
Fig. 5.3: Snapshots of the primary system showing spatial distribution of ions at (a) t =
0.1 ps before the shock front strikes the nanobubble, (b) t =1.7 ps when the water
hammer shock forms, and (c) t = 3.5 ps when the number of ions has reached a
stationary value. Positive and negative ions are shown in red and blue, respectively.
................................................................................................................................... 81
Fig. 5.4: Ions and carbon atoms in the primary system at t = 3.5 ps. (a) Side view. (b)
Front view. Positive and negative ions are shown in red and blue, respectively.
Green spheres depict carbon atoms........................................................................... 82
Fig. 5.5: Four dominant pathways to splitting of H
2
O molecules observed in MD
simulations. Each row displays one pathway, where the reactants and intermediate
states shown respectively in the left and right panels. Carbon, hydrogen and oxygen
atoms are shown respectively in green, white and red.............................................. 83
Fig. 5.6: A pathway for formation of a formic acid molecule. Each row displays a
reaction, left and right panels being reactants and products, respectively. These two
reactions take place after the reaction shown in Fig. 6.5(c). The formic acid
molecule is seen in panel (d). Transfers in panel (a) are: O (2) and H (3) to C (1), H
(4) to O (5), H (6) to O (7), H (8) to O (9). Moreover, unstable bond between H (10)
and H (11) breaks. In panel (c), H (2) transfers to O (1) and bond between O (1) and
C (3) breaks............................................................................................................... 84
Fig. 5.7: (a) Number of CO2 molecules as a function of time in the primary system. At t
= 1.4 ps, H2O-CO2 reactions start, and some unstable products start decaying at t =
2.3 ps. (b) Front view (the direction of shock propagation is perpendicular to the
figure and is outward) of carbon atoms and oxygen molecules at the end of
simulation at t = 3.5 ps in the primary system. Green spheres depict carbon atoms,
and oxygen molecules are shown in red. .................................................................. 86
Fig. 5.8: Formation of an oxygen molecule from an H
2
O
2
molecule. Carbon, hydrogen
and oxygen atoms are shown in green, white, and red, respectively. Narrow sticks
between oxygen and hydrogen atoms indicate hydrogen bonds............................... 87
Fig. 5.9: Molecules involved in formation of an oxygen molecule. Carbon, hydrogen, and
oxygen atoms are shown in green, white, and red, respectively. Narrow bars between
oxygen and hydrogen atoms indicate hydrogen bonds. Here, a 2-carbon molecule is
surrounded by 2 water molecules and a water dimer. Oxygen atoms marked by 1 and
2 separate from their molecules and form an oxygen molecule. .............................. 88
ix
Abstract
The shock-induced collapse of nanobubbles in water is investigated using
molecular dynamics simulations based on a reactive force field. Monitoring the collapse
of a cavitation nanobubble, we observe a focused nanojet at the onset of bubble shrinkage
and a water hammer shock wave upon bubble collapse. The nanojet length scales linearly
with the nanobubble radius, as observed in experiments on micron-to-millimeter size
bubbles. The shock induces dramatic structural changes, including an ice-VII-like
structural motif at a particle velocity of approximately 1 km/s. The incipient ice VII
formation and the calculated Hugoniot curve are in good agreement with experimental
results. Moreover, a substantial number of positive and negative ions appear when the
nanojet hits the distal side of the nanobubble and the water hammer shock forms.
Furthermore, two promising applications of shock-induced nanobubble collapse have
been explored. Our simulations of poration in lipid bilayers due to shock-induced
collapse of nanobubbles reveal penetration of nanojets into lipid bilayers. The nanojet
impact generates shear flow of water on bilayer leaflets and pressure gradients across
them, which transiently enhance the bilayer permeability by creating nanopores through
which water molecules translocate across the bilayer. The effects of nanobubble size and
temperature on the porosity of lipid bilayers are examined. Finally, the shock-induced
collapse of CO
2
-filled nanobubbles in water is investigated. The energetic nanojet and
high-pressure water hammer shock formed during and after collapse of the nanobubble
trigger mechano-chemical H
2
O-CO
2
reactions, some of which lead to splitting of water
x
molecules. The dominant pathways through which splitting of water molecules occur are
identified.
1
1 Introduction
1.1 Outlook
The research presented in this dissertation consists of simulations of three related
systems. Before explaining these systems and their connection, let us first pause on the
nature of these simulations. Computer simulations are performed mostly for two
purposes: 1) testing a theory or model and 2) invoking a theory or model to predict some
properties of a system.
A computer simulation can provide a validation for a theory. For example, the
agreement between the experimental mass of the proton and the simulated value served
as a test for the new lattice quantum field theory. Similarly, a computer simulation can be
used to examine a model. As an example, an astrophysicist can perform a simulation to
check a model of the sun’s outer atmosphere against the data experimentally available to
find how well the model works.
Alternatively, a computer simulation can be utilized for predictions based on a
theory or model. Having checked his model against available data, our astrophysicist can
perform new simulations to predict the next coronal mass ejection. As another example,
properties of materials can be predicted especially in extreme regimes in which
experiments are hard or impossible to perform provided a suitable model is available.
2
Although most computer simulations are performed for the two above mentioned
purposes, they also can be used for purely exploratory purposes. Clearly, discoveries in
its fundamental sense are not possible, since the outcomes of a computer simulation
depend on the input theory or model. However, one can explore the possibilities of some
scenarios or phenomena via computer simulations. For example, in 1950s, statistical
physicists were wondering whether spherical particles with stiff short-range repulsion and
without any attraction can form a crystal or not. In a famous simulation work, Alder and
Wainwright (Alder and Wainwright 1957) and Wood and Jacobson (Wood and Jacobson
1957) showed such a system has a first-order freezing transition. Simulations presented in
this dissertation are also exploratory simulations. The next two sections discuss which
possibilities are explored and why they are important.
1.2 Shock-induced Collapse of Nanobubbles
Cavitation bubbles form naturally and gas bubbles can be synthesized. Instability
of bubbles has been studied analytically (Plesset and Chapman 1971; Prosperetti 1977).
Of more interest are analytical (Hilgenfeldt, Lohse et al. 1996; Hao and Prosperetti 1999;
Lin, Storey et al. 2002) and experimental (Holt and Gaitan 1996; Gaitan and Holt 1999;
Ketterling and Apfel 2000; Ketterling and Apfel 2000) studies of violent collapse of
bubbles. Experimental studies on microbubbles reveal that the impact of a pulsed shock
on the proximal side of the bubble shrinks and accelerates it in the direction of the shock
propagation. Furthermore, the liquid around the collapsing bubble forms a jet, which
creates a protrusion and water hammer shock wave when it hits the distal side of the
3
bubble. Later, the jet breaks up and in addition to the bubble, a new tiny bubble also
forms; see Fig. 1.1.
Fig. 1.1: Three stages of shock-induced collapse of a microbubble are shown
schematically in the left panel. When the shock front strikes on the bubble, it shrinks
because of outside overpressure, and center of the bubble accelerates upward. At this
time, liquid surrounding the bubble enter into it and forms a jet, which is focused towards
the distal point of the bubble. Finally, the jet breaks and leaves behind a tiny bubble.
Photographs of the bubble before and after passage of the shock wave are shown in the
right panel. Picture courtesy of (Ohl and Ikink 2003).
The behavior of a collapsing bubble depends on its location relative to objects
embedded in the fluid as well as the shock wave characteristics (amplitude and pulse
width) and the initial size of the bubble. Shock-induced bubble collapse phenomena have
been studied in the context of a single bubble (Brujan, Keen et al. 2002; Tomita and
Kodama 2003; Sankin and Zhong 2006) as well as multiple bubbles (Tufaile and
Sartorelli 2002) near a rigid boundary (Tomita and Shima 1986; Philipp and Lauterborn
1998), elastic membrane (Yamanoi and Tamagawa 2006; Ohl, Klaseboer et al. 2009), and
in confined environments (e.g., microfluidic or lab-on-a-chip systems) (Ajaev and Homsy
2006; Zwaan, Le Gac et al. 2007).
4
Jetting and water hammer shocks can cause significant damage in materials. This
problem is encountered in disintegration of blades of ship propellers, pipelines and pump
blades (Brennen 1995). Moreover, experimental (Kodama and Takayama 1998) studies
and simulation works (Johnsen and Colonius 2008) shows that bubble collapse can cause
damage to the kidney tissue in shock wave extracorporeal lithotripsy.
Motivated by aforementioned studies on bubble collapse and lack of any
experiment on collapse of nanobubbles, we studied shock-induced collapse of
nanobubbles in water. The goals were to study the structure of water under shock in
atomistic level and to study the dynamics of the collapse. More specifically, we were
interested in discovering whether or not jetting phenomenon can happen in nanometer-
sized bubbles and to what extent nanobubble collapse is violent. Answers to these
questions will be discussed in detail in Chapter 3.
An interesting aspect that neither experiments nor simulations have yet examined
is the mechano-chemistry of collapse of bubbles. It is well known that water molecules in
the liquid phase dissociate into hydronium (H
3
O
+
) and hydroxide (OH
-
) ions. These
autoionization events in liquid water are rare, occurring once in approximately 10 hours.
This, naturally, raises the question: Does the collapse of a cavitation bubble enhance
autoionization in liquid water? The results presented in Chapter 5 answer this question.
1.3 Applications of Shock-induced Collapse of Nanobubbles
In recent years, non-invasive drug- and gene-delivery approaches have garnered
significant interest because of direct applications in cancer treatment and gene therapy.
5
Much of the experimental effort is focused on designing a targeted approach that has both
spatial and temporal specificities. Research in this area relies mainly on the use of electric
fields or pressure waves to enhance the permeability of cell membranes. In one of the
most commonly used techniques known as electroporation, electric fields are applied
across the cell to increase the cell membrane permeability (Vernier, Sun et al. 2006;
Golberg and Rubinsky 2010). Reversible electroporation, in which the cell permeability
is enhanced temporarily, is used for drug delivery and gene therapy. Electric fields
applied over a sufficiently long time can kill the cell because of temperature elevation
resulting from Joule heating. This irreversible electroporation process is commonly used
in the food industry to inactivate microbes and also in minimally invasive treatment of
cancerous tissues.
Sonoporation is another promising DNA-, protein- and drug-delivery approach.
To achieve high efficiency in sonoporation, in vivo gas bubbles are used in conjunction
with diagnostic level ultrasound exposures (Rosenthal, Sostaric et al. 2004; Mitragotri
2005; Yamanoi and Tamagawa 2006; Newman and Bettinger 2007; Zhou, Cui et al.
2008; Kudo, Okada et al. 2009). In sonoporation approach, the collapse of microbubbles
generates liquid jets and radial spreading flows (Ohl, Arora et al. 2006) that can make the
cell membrane transiently permeable to molecular entry. While the aforementioned
experimental studies investigate the effect of microbubbles in sonoporation, we have
performed simulations to study effects of shock-induced collapse of nanobubbles in
permeability of cell membranes. We showed that nanobubbles are big enough to induce
6
transient nanopores in cell membranes. Moreover, we investigated effects of nanobubble
size and temperature, which will be discussed in chapter 4.
The sun is by far the largest resource of energy available to us. For years,
photovoltaic cells have been in use for the generation of electricity from solar radiation.
Improving efficiency and utilizing cheaper materials resulting in lower cost of generated
power are continuing challenges for researchers in this area (Lewis and Nocera 2006). On
a large scale, relatively recent concentrating solar power plants focus solar radiation on
small target areas and generate heat at temperatures ranging from 200-1000ºC, which in
turn can be converted to electricity through conventional steam cycles (Roeb and Muller-
Steinhagen 2010). The storage problem arising from the intermittent nature of solar
energy is the major challenge here. Pumped hydropower and molten salt are used to store
electricity, and new materials such as special concretes or high-temperature phase-change
materials are emerging.
Chemical fuels are one of the possible solutions for the storage problem.
Moreover, they can substitute fossil fuels, commonly used in transportation. The desire
for sustainable production of chemical fuel sufficient for the current global demand can
be realized by production of hydrogen from H
2
O and CO
2
as feedstock, which raises the
challenge of efficient H
2
O and CO
2
splitting (Turner 2004). Direct thermolysis of H
2
O
and CO
2
requires operating temperatures above 2000 °C. At present, two-step
thermochemical cycles using metal oxides are promising approaches to reduce the
required operating temperature (Steinfeld 2005; Kodama and Gokon 2007). In this
approach, the metal oxide reduces in the first step, and then oxidizes in the second step
7
establishing a cycle. Different metal oxides and reactor designs have been examined.
Low conversion efficiencies and usage of precious materials are common challenges.
Ferrite-based oxides have high reduction temperatures and show relatively slow reaction
rates—which becomes smaller by time because of sintering (Perkins and Weimer 2009).
Volatile oxides require rapid quenching of gaseous products to avoid recombination
(Gstoehl, Brambilla et al. 2008; Miller, Allendorf et al. 2008; Schunk, Lipinski et al.
2009; Schunk and Steinfeld 2009). Cerium oxide (ceria) is a new promising oxide, which
despite its non-stoichiometric redox cycle exhibits a high fuel production rate (Kaneko,
Miura et al. 2007; Miller, Allendorf et al. 2008; Chueh and Haile 2009; Abanades, Legal
et al. 2010; Chueh and Haile 2010). Furthermore, alternative reactor designs exposing
ceria directly to concentrated solar radiation has been proposed (Chueh, Falter et al.
2010).
Recently, an alternative approach for the direct conversion of mechanical energy
to chemical energy has been proposed, which utilizes vibrating piezoelectric microfibers
to decompose water molecules (Hong, Xu et al. 2010). On the other hand, hot
temperature and pressures achieved during the collapse of cavitation bubbles are known
to induce chemical reactions (Didenko and Suslick 2002; Kuijpers, van Eck et al. 2002;
Flannigan and Suslick 2005). Similarly, as mentioned, collapses of nanobubbles are
accompanied by energetic nanojets and high-pressure water hammer shock waves.
Putting all together, a natural question to ask is whether shock-induced collapse of a
nanobubble can be utilized to split H
2
O and CO
2
molecules. Chapter 5 presents results of
the simulations performed to answer this question.
8
2 Reactive Force Field Molecular Dynamics
Simulations on Parallel Computers
2.1 Molecular Dynamics Simulations
Molecular dynamics (MD) can be described as brute force integration of
Newton’s equations of motion to obtain the phase space trajectories of classical many-
body systems. Here, “classical” implies that each particle in the system experiences a
force, which determines its motion according to Newton’s second law. However, one can
also construct effective forces that are designed to capture some quantum mechanical
properties of the system. For example, there are reactive force fields, which describe
chemical bond formation and dissociation (Johnston and Parr 1963; Root, Landis et al.
1993; van Duin, Dasgupta et al. 2001). Apparently, forces constructed in this way are
system dependent and should be tuned (or even reconstructed) when considering new
systems.
2.1.1 Time Integration
There are numerous time-integration algorithms devised for MD simulations,
among which the velocity Verlet algorithm is widely used. Tuckerman’s (Tuckerman,
Berne et al. 1992) derivation of the velocity Verlet algorithm from the Liouville
formulation of classical mechanics shows that this algorithm is time-reversible and
9
preserves phase-space volume. An outline of Tuckerman’s approach (Frenkel and Smit
2002) is given in the following.
Let us consider a function , which depends on the momenta and
coordinates of N particles in a classical many-body system, but does not depend on time
explicitly. Here, 3N-dimensional vectors p and r denote momenta and coordinates of all
particles, collectively. The time derivative of f is
€
˙
f = ˙ r
∂f
∂r
+ ˙ p
∂f
∂p
≡ iLf ,
(2.1)
where the Liouville operator L is defined as
€
iL = ˙ r
∂
∂r
+ ˙ p
∂
∂p
. The above equation can be
integrated formally to yield
(2.2)
Let us also consider the following Trotter identity:
€
e
(A +B)
=
P→∞
lim
(e
A /2P
e
B /P
e
A /2P
)
P
= (e
A /2P
e
B /P
e
A /2P
)
P
e
O(1 P
2
)
,
(2.3)
where A and B are two operators that do not necessarily commute. In MD simulations,
is known and the goal is to find using a number of MD steps
each of length Δt. Identifying and
(2.4)
(2.5)
10
each single MD time step is equivalent to applying the following operator
(2.6)
which consists of three elementary operators. Let us consider the effect of each
elementary operator on coordinates and momenta separately. Starting from the right, the
application of the operator shifts the momenta and leaves the coordinates
unchanged.
(2.7)
Similarly, the application of the middle elementary operator to the result of the
previous operation shifts the coordinates and leaves the momenta unaffected.
(2.8)
Finally, application of the left elementary operator on the result of the
previous operation yields
(2.9)
The overall effect of the above three elementary transformations is equivalent to the
velocity Verlet algorithm; see Table 2.1, where r, p, and F denote coordinate, momentum,
and force vectors for each individual particle.
11
Table 2.1: Velocity Verlet Algorithm
Calculate F and update p as
Update r as
Calculate F again and update p as
As mentioned above, the velocity Verlet algorithm has the advantage of being
time-reversible and phase-space-volume preserving. The former is desired, since classical
mechanics formulation is time-reversible, and the latter is important, since it prevents
from a drift in energy, and results in a more stable algorithm. Time-reversibility follows
from the fact that past and future coordinates and momenta appear symmetrically in the
algorithm. Now, let us consider the Jacobian of the transformation from
to . The Jacobian of the overall transformation equals the product of
Jacobians of the three elementary transformations. Equations (2.7-9) show that in each
elementary transformation, only one of p or r shifts. Moreover, shifts in p depend on only
r (
€
˙ p =F(r)) and the shift in r depends on only p (
€
˙ r = p). Hence, Jacobians of the three
elementary transformations and consequently the overall transformation are unity. The
overall Jacobian being unity implies that the velocity Verlet algorithm preserves phase-
space-volume.
2.1.2 Periodic Boundary Conditions and Minimum Image Convention
Each MD box has boundaries, and consequently surface effects are present. To
simulate properties of bulk materials, surface effects should be eliminated. Imposing
12
periodic boundary conditions (PBC) is an approximation method to simulate bulk
materials. In this approximation, replications of the MD box fill the entire space thus
effectively, forming a bulk material. Fig. 2.1 schematically depicts a 2-dimensional MD
box (gray cell) and its nearest-neighbor images. Since all cells are the same, it is
sufficient to keep the coordinates of particles only in the MD box. Whenever a particle
moves, all its images move the same way; therefore, if a particle exits the MD box (blue
particle), its image particle (red particle) enters the MD box. Evidently, approximating a
bulk material by replication of a small system (the MD box) introduces some issues,
which depend both on the phenomenon under consideration and the range of inter-
particle potentials (Allen and Tildesley 1987). However, based on common experience in
MD simulations, if the potential range is such that a particle and its images interact
negligibly, and if we limit ourselves to regimes away from phase transitions, there is little
effect on thermodynamic and structural properties arising from imposition of PBC.
13
Fig. 2.1: Schematic illustration of periodic boundary condition in two dimensions. The
MD cell is marked by gray, and its nearest-neighbor images are shown. If a particle exits
the MD box (blue particle), its image particle (red particle) enters the MD box
Introducing PBC eliminates surface effects, but brings a practical problem. With
PBC, the number of particles is effectively infinite, and thus force calculations become
impractical. Metropolis et al. (Metropolis, Rosenbluth et al. 1953) introduced another
approximation called “minimum image convention” to solve this problem. Minimum
image convention hypothesizes that each particle interacts with another particle if they
are not apart more than half of the length of the MD box; see Fig. 2.2. Let us consider a
particle i (shown in red), and partition all particles in classes such that all particles in a
class are images of a single particle in the MD box. One of the classes is the class of
particles, which are images of particle j. Among all particles in this class (shown in green
and blue), according to MIC, only the one being the nearest to particle i (shown in blue)
will interact with particle i.
14
Fig. 2.2: Schematic illustrationof the minimum image convention in two dimensions.
Particle i does not interact with particles shown in green, and only interacts with the
particle shown in blue.
2.1.3 Potential Cutoff Distance and Linked-list Cell Method
Let us consider a MD box containing N particles and let us focus on pairwise
interactions with the minimum image convention. The force on each particle comprises
contributions from all other particles. Therefore, force calculation for each particle
requires O(N) operation. Considering that force calculation is the most expensive part of
a MD program, the MD simulation has a complexity of O(N
2
). Introduction of a potential
cutoff distance is another approximation approach for reducing the complexity of MD
simulations for short-range potentials. Short-range potentials like the Lennard-Jones
potential fall rapidly as the distance between a pair increases; hence, by choosing a good
cutoff distance (r
c
),
neglected contributions do not result in a significant error.
15
Mere introduction of the cutoff distance does not reduce the complexity of the
MD program, since still distances between all pairs should be calculated and compared
with the cutoff distance. However, cutoff distance used together with the cell index
method (Quentrec and Brot 1973; Hockney and Eastwood 1981) or its memory-efficient
version called “inked-list cell method” (Knuth 1973; Hockney and Eastwood 1981)
reduces the complexity of the MD program to O(N). Invoking the linked-list cell method,
the MD box is divided into N
c
cells with dimensions of , which are slightly larger than
the cutoff distance. Here, α denotes Cartesian coordinate directions x, y, and z. In this
method, a linked-list data structure keeps track of particles in each cell as follows. First,
all cells are given a sequential index. Then, two one-dimensional arrays are created as
shown in Fig. 2.3. The first array head is of size N
c
and the second array lcsl is of size N.
Table 2.2. shows the way these two arrays are constructed. Let us consider how particles
inside a cell can be retrieved through an example using Fig. 2.3. To retrieve particles in
cell 3, the first atom is found in head[3], which is 7; the second atom is found in lcsl[7]
which is 3; the third atom is found in lcsl[3] which is 1; finally, lcsl[1] is empty implying
that all atoms in cell 3 are retrieved.
When calculating the force on a particle i in a cell with index k, by utilizing
linked-list cell method, only particles inside the cell k and neighboring cells are
considered. All other particles have a distance larger than r
c
as measured from the
location of the particle i. In three dimensions, each cell has 26 neighbors, and number of
particles in each cell is O(1). Accordingly, the complexity of force calculation reduces to
O(N) as compared to O(N
2
) for brute force approach.
16
Table 2.2: Pseudocode for constructing head and lscl arrays
and
for each atom i
find the cell k to which atom i belongs
end for
Fig. 2.3: Data structure for linked list method. Here, “E” stands for Empty. Picture
courtesy of Aiichiro Nakano.
2.1.4 Shifted and Shifted-Force Potential for Pairwise Interactions
Introducing a cutoff distance r
c
for pairwise interactions expedites potential and
force calculation, but sets forth a discontinuity in potential at the cutoff distance. The first
attempt to solve this problem was to define a shifted potential as
17
(2.10)
where r
ij
is the distance between particles i and j, and V
c
denotes the potential at the cutoff
distance. Shifted potential resolve discontinuity in potential, so force can be defined.
However, the force will not be continuous resulting in instability in numerical solutions
of the equations of motions (Stoddard and Ford 1973; Nicolas, Gubbins et al. 1979;
Powles, Evans et al. 1982). The shifted-force potential defined as
(2.11)
removes the discontinuity in force by adding a small linear term to the potential. The
shifted-force potential goes smoothly to zero at the cutoff distance, and therefore
instability problems do not rise.
2.1.5 Molecular Dynamics in Various Ensembles
So far, we have talked about MD in a microcanonical ensemble i.e., constant
number of particles, constant volume, and constant energy (since Newton’s equations of
motion preserve total energy in the absence of external forces). NVE ensemble is actually
what is referred to as microcanonical ensemble in statistical mechanics. Moreover, we
showed that the velocity-Verlet algorithm preserves phase-space volume. To put it
another way, if we consider a system at time t at a phase space point Γ, when this system
leaves Γ at time t+1 (let time be discrete), another system arrives at Γ. Therefore, the
motion of systems in phase space resembles a stream. Actually, there might be several
18
streams passing through different regions of the phase space. However, if we assume that
there is only one stream passing through all points in phase space (i.e., if the system is
ergodic), then we can conclude that any system initially at a specific phase space point
will eventually pass through all points accessible to that system (points on a constant
energy shell). In this case, the time average and ensemble average will be the same.
When performing an MD simulation in the microcanonical ensemble, we start
with an arbitrary state in microcanonical ensemble and assuming a phase-space-
preserving time integration algorithm and also assuming ergodicity, we sample the
microcanonical ensemble by following the trajectory of the system obtained from
integration of Newton’s equation of motion. In practice, our simulation time is limited, so
there is no guarantee that we have sampled all states in the ensemble. We usually
continue the simulation until the thermodynamic quantity of interest converges up to a
desired precision.
In statistical mechanics, different ensembles such as canonical and grand
canonical ensembles are considered. Since, ensembles are basically artificial constructs,
one expects to obtain the same value for observables in thermodynamic limit. This
expectation is believed to be true for ensembles commonly used in statistical mechanics
provided we are not close to phase transitions (Fisher 1964). There are two widely used
methods for calculating the value of an observable in an ensemble other than
microcanonical (e.g., canonical): Monte Carlo (stochastic) and MD (deterministic)
simulations. In the simulations presented in this dissertation, we do not use the Monte
19
Carlo technique, but focus on MD simulations. When thinking about non-NVE MD
simulations, two questions arise:
1) How can one perform an MD simulation in ensembles other than the
microcanonical ensemble?
2) Considering finite numbers of particles in any simulation (as opposed to infinite
number of particles in the thermodynamic limit), are values of observables
obtained from simulations in different ensembles the same?
The answer to the second question can be found in (Lebowitz, Percus et al. 1967). Let us
focus on the first question. To be able to perform an MD simulation in an ensemble other
than microcanonical, we need a prescription of how to move a system initially at a
specific point of the phase space to other points accessible for the system in the ensemble
under consideration. To be useful, this prescription should satisfy some conditions:
- The prescription should be phase-space-preserving for the same reason discussed
above for microcanonical ensemble.
- We should be able to argue that ergodicity holds, even though we cannot prove it.
- Since simulation time is limited and since different states in the ensemble under
consideration (which is not a microcanonical ensemble) have different
probabilities, we should pass through more probable states sooner. In other words,
the average of desired observable should converge during the simulation.
In the following, we briefly mention the prescriptions to perform MD simulations in a
canonical ensemble (i.e., constant number of particles, volume, and temperature), which
has been used in simulation to be presented later.
In 1984, Nosé showed that one can perform deterministic MD at constant
temperature (Nose 1984; Nose 1984). His approach is based on an extended Lagrangian
20
with artificial coordinates and velocities. Later, it became common to use the Nosé
scheme in the formulation of Hoover, which is referred to as the Nosé-Hoover thermostat
(Hoover 1985; Hoover 1986). Then, it was noticed that there are exceptional cases (e.g., a
simple harmonic oscillator) for which the Nosé-Hoover thermostat does not generate a
canonical distribution. In 1992 Martyna et al. solved the problem by coupling a number
of thermostats together (Martyna, Klein et al. 1992). The equations of motion for a
system consisting of N particles with M Nosé-Hoover chains are given by
€
˙ r
i
=
p
i
m
i
, (2.12)
€
˙ p
i
=−∇
r
i
U r
i
{ }
( )
− p
i
p
s
1
Q
1
, (2.13)
€
˙ s
k
=
p
s
k
Q
k
for k =1,...,M, (2.14)
€
˙ p
s
1
=
p
i
2
m
i i=1
N
∑
− 3Nk
B
T
eq
⎡
⎣
⎢
⎤
⎦
⎥ −
p
s
2
Q
2
p
s
1
, (2.15)
€
˙ p
s
j
=
p
s
j−1
2
Q
j−1
− k
B
T
eq
⎡
⎣
⎢
⎢
⎤
⎦
⎥
⎥
−
p
s
j+1
Q
j +1
p
s
j
for j = 2,…,M −1 , (2.16)
€
˙ p
s
M
=
p
s
M−1
2
Q
M−1
− k
B
T
eq
, (2.17)
where and are thermostat coordinates and momenta, respectively. k
B,
and T
eq
are the Boltzmann constant and the desired temperature, respectively. Q
k
is the effective
mass of the kth thermostat.
21
2.1.6 Physical Quantities
According to equivalence of ensembles, basic thermodynamic quantities can be
calculated as averages in any convenient ensemble provided a suitable phase function
exists. Kinetic, potential, and total energies can be calculated as
(2.18)
where , , and are relevant phase functions defined by
€
V = V
1
i
∑
(r
i
) + V
2
(r
i
j>i
∑
,r
j
) + V
3
k> j>i
∑
j>i
∑
i
∑
i
∑
(r
i
,r
j
,r
k
)…, (2.19)
(2.20)
(2.21)
Here, r
i
and p
i
denote the coordinate and momentum vector of particle i, respectively. V
1
,
V
2
, and V
3
are respectively 1-body, 2-body, and 3-body potential energy functions.
Temperature and pressure can be calculated using the virial theorem written in the
form of independent generalized coordinates q
k
and generalized momenta p
k
as
(2.22)
where k
B
is the Boltzmann constant. It is common to define an instantaneous “kinetic
temperature” whose average equals thermodynamic temperature T. If momenta
appear as squared terms in the Hamiltonian of a system, and if there are N
c
constraints
(e.g. molecular constraints or constraints on the motion of the center of mass),
€
T =
1
(3N−N
c
)k
B
p
i
2
m
i i
∑
. (2.23)
22
A system under shock is not in thermal equilibrium and thermodynamic temperature is
not meaningful for such a system. However, there is a way to define “temperature” in
shock simulations (Holian 1995).
Considering Hamilton’s equation
€
∂H
∂q
k
=− ˙ p
k
, (2.24)
together with the virial theorem, we arrive at the following relation
€
1
3
r
i
⋅ f
i
tot
i
∑
=
1
3
r
i
⋅ f
i
ext
i
∑
+
1
3
r
i
⋅ f
i
int
i
∑
=−Nk
B
T, (2.25)
where , , and are internal, external, and total force vectors on particle i. The
contribution from the external component of the force can be identified as
(2.26)
where P and V are the thermodynamic pressure and volume of the system, respectively.
By rearranging,
(2.27)
where ρ is number density. Defining an instantaneous pressure (with the average of
P),
€
P =
1
3V
r
i
⋅ f
i
int
i
∑
+ρk
B
T. (2.28)
In simulations with PBC, it is essential to note that the above equation is not well defined.
For pairwise interactions, invoking to Newton’s third law,
23
(2.29)
We should note that defining local pressure based on the virial theorem is controversial.
To define local pressure with the virial theorem, the sum over atoms in Eq. (2.28) should
be broken into contributions from individual atoms, which makes V ambiguous. However,
there is a way to by pass this subtlety and define local pressure in MD simulation more
rigorously (Tsai 1979; Cheung and Yip 1991; Sun, Wang et al. 2006).
Finally, let us consider a structural quantity: the radial distribution function g(r).
Formally, g(r) can be defined by integrating configurational distribution function over all
coordinates except two. Incorporating a normalization factor (McQuarrie 1976; Hansen
and McDonald 1986), in the canonical ensemble,
€
g(r) =
N(N −1)
ρ
2
Z
NVT
d
∫
r
3
dr
4
dr
N
exp(−βV (r
1
,r
2
,r
N
)), (2.30)
where and N, ρ, and Z
NVT
are number of particles, number density, and
partition function in the canonical ensemble, respectively. In practice, an equivalent
definition is used in MD simulations,
(2.31)
where . To implement the above formula, a histogram with bins of length δr
is considered, which keeps the number of pairs with separation between and
in the bin (i+1).
24
2.2 ReaxFF: A Reactive Force Field
Shock-induced nanobubble collapse simulations (with the exception of
simulations performed for poration of lipid bilayers) utilized a reactive force field,
ReaxFF (van Duin, Dasgupta et al. 2001; Strachan, van Duin et al. 2003; van Duin,
Strachan et al. 2003). In ReaxFF, energy of a system depends on atomic coordinates,
bond order (BO) between atomic pairs, and atomic charges.
2.2.1 Bond Order Calculation
In ReaxFF, bond order BO
ij
between a pair of atoms i and j describes the nature of
covalent bonds between the pair and is calculated in two steps. In the first step, raw bond
orders BO
ʹ′
ʹ′
ij
are computed as a function of interatomic distances,
(2.32)
where is the interatomic distance between atoms i and j (r
i
denotes the
position of atom i). In the second step, raw bond orders are corrected to yield better
atomic valences as
(2. 33)
(2.34)
(2.35)
25
(2.36)
(2.37)
(2.38)
(2.39)
€
f
3
(
ʹ′
Δ
i
,
ʹ′
Δ
j
) =−
1
p
boc2
ln
1
2
exp −p
boc2
ʹ′
Δ
i
( ) +exp −p
boc2
ʹ′
Δ
j ( ) [ ]
⎧
⎨
⎩
⎫
⎬
⎭
, (2.40)
(2.41)
(2.42)
(2.43)
where n(i) is the set of neighbor atoms within single covalent bond cutoff distance for
atom i. All parameters in equations of this section are listed in Table 2.3. When bond
orders are corrected, corrected overcoordination of atom i is computed as
(2.44)
26
Table 2.3: Parameters in bond order functions
bond parameters
bond radius parameters
valencies of atom i
overcoordination parameters
1-3 bond order corrections
2.2.2 Charge Calculation
In ReaxFF, atomic charges are obtained using an electronegativity equalization
method (EEM) (Mortier, Ghosh et al. 1986; Janssens, Baekelandt et al. 1995).
Sanderson’s electronegativity equalization principle (Sanderson 1976; Parr and Pearson
1983; Sanderson 1983) postulates that all constituent atoms in a molecule should have
equal electronegativity. Moreover, density functional theory shows (Mortier,
Vangenechten et al. 1985; Mortier, Ghosh et al. 1986; Vangenechten, Mortier et al. 1987)
that electronegativity of atom i in a molecule can be calculated by
(2.45)
€
A
i
=χ
i
0
+Δχ
i
,
(2.46)
(2.47)
where N is the number of atoms in the molecule.
€
q
i
, , and are charge,
electronegativity, and hardness of atom i, respectively.
€
R
ij
denotes the distance between
atoms i and j. κ, , and are empirical parameters, which describe the molecular
27
environment. According to Sanderson’s electronegativity equalization principle,
electronegativity of all atoms in a molecule is equal to molecule’s electronegativity,
(2.48)
where is molecular electronegativity. Conservation of charge yields
(2.49)
where Q is the total charge of the molecule. Finally, atomic charges can be obtained in
terms of , , κ, and Q, as
€
B
1
κ
R
12
κ
R
1N
−1
κ
R
21
B
2
κ
R
2N
−1
κ
R
N1
κ
R
N2
B
N
−1
1 1 1 0
⎡
⎣
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎤
⎦
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
q
1
q
1
q
1
χ
⎡
⎣
⎢
⎢
⎢
⎢
⎢
⎢
⎤
⎦
⎥
⎥
⎥
⎥
⎥
⎥
=
−A
1
−A
2
−A
N
Q
⎡
⎣
⎢
⎢
⎢
⎢
⎢
⎢
⎤
⎦
⎥
⎥
⎥
⎥
⎥
⎥
.
(2.50)
In our implementation of ReaxFF, we minimize the electrostatic energy with a
charge neutrality constraint, which is equivalent to EEM method. Let energy be an
arbitrary function of atomic charges, . The charge neutrality constraint reads
. For minimizing energy subject to the charge neutrality constraint, we can use
Lagrange multiplier method. Therefore, can be minimized without
a constraint where is the Lagrange multiplier. When is minimized,
(2.51)
28
Noting that is the chemical potential of atom i, or negative of electronegativity of
atom i, then we can see that all atoms have the same electronegativity. We use conjugate
gradient method to minimize energy.
2.2.3 Potential Energy
In ReaxFF, the potential energy of a system consists of twelve components:
(2.52)
where bonding energy E
bond
, lone pair energy E
lp
, overcoordination energy E
over
,
undercoordination energy E
under
, valence-angle energy E
val
, penalty energy E
pen
, 3-body
conjugation energy E
coa
, torsion-angle energy E
tors
, four-body conjugation energy E
conj
,
and hydrogen bonding E
hbond
are energy of valence interactions while van der Waals
energy E
vdWaals
and Coulomb energy E
Coulomb
are energy of non-covalent interactions. The
interactions in the above equation can be grouped into 1-, 2-, 3-, and 4-body interactions
according to number of explicitly involved atoms; see Table 2.4 and Fig. 2.4. In ReaxFF,
energies of valence interactions depend on bond orders, which ensure smooth vanishing
of these energies as bonds break.
29
Table 2.4: Classification of interactions in ReaxFF
Number of explicitly involved atoms
atoms
Energy terms
1 E
lp
,
E
over
, E
under
2 E
bond
, E
vdWaals
, E
Coulomb
3 E
val
,
E
pen
,
E
coa
,
E
hbond
4 E
tors
, E
conj
Fig. 2.4: Schematic illustration of atomic configurations in energy terms: (a) 1-body, (b)
2-body, (c) 3-body, (d) 4-body, (e) hydrogen-bonding and (f) non-covalent iterations,
respectively. A red sphere represents the position of the primary atom in each energy
term. Open bars represent covalent bonds, while dotted lines are non-covalent bonds.
Picture courtesy of (Nomura, Kalia et al. 2008).
2.2.3.1 1-body Interactions
Lone pair energy, over coordination energy, and undercoordination energy are
defined for single atoms. Number of lone pairs around atom i is defined as
€
n
lp,i
=
Δ
i
e
2
⎢
⎣
⎢
⎥
⎦
⎥
+exp −p
lp1
2 +Δ
i
e
−2
Δ
i
e
2
⎢
⎣
⎢
⎥
⎦
⎥
⎛
⎝
⎜
⎞
⎠
⎟
2
⎢
⎣
⎢
⎢
⎥
⎦
⎥
⎥
, (2.53)
30
where measures the difference between the total number of outer shell electrons (e.g. 6
for oxygen, 4 for silicon, or 1 for hydrogen) and the sum of bond orders around atom i,
(2.54)
As the total bond order associated with atom i starts to exceed its normal value, Eq.
(2.53) allows a lone pair to gradually break up causing a deviation from the optimal
number of loan pairs
€
n
lp−opt,i
defined as
. (2.55)
Once is known, the penalty energy function for atoms with lone electron pairs
exceeding the optimal number is given by
(2.56)
Even after correction of bond orders as described in §2.2.1, some atoms may still
have bond orders exceeding what the theory of valence electrons predicts. The remaining
overcoordinations receive an energy penalty. However, as mentioned before, Eq. (2.53),
allows lone pair electrons to break up and the resulting overcoordinations should not be
punished. Therefore, first we calculate corrected overcoordination as
(2.57)
After obtaining corrected overcoordinations, each over coordinated atom i will receive an
energy penalty given by
31
(2.58)
If an undercoordinated atom i is connected to another undercoordinated atom j
and if bonds between i and j have partly π-bond character, the energy contribution arising
from resonance of π-electrons between attached undercoordinated atomic centers will be
taken into account as
(2.59)
Parameters appearing in 1-body energy functions are listed in Table 2.5.
Table 2.5: Parameters in 1-body energy functions
valency of atom i
normal number of lone pair electrons for atom i
parameters
lone pair parameter
lone pair energy parameter
overcoordination parameters
undercoordination energy parameter
undercoordination parameters
σ-bond energy parameter
32
2.2.3.2 2-body Interactions
Bond energy, van der Waals energy, and Coulomb energy are defined for all pairs
of atoms provided the distance between the pair is less than cutoff length of the
corresponding interaction. Bond energy E
bond
sums energy contributions from σ-, π-, and
ππ-bonds,
(2.60)
As mentioned above, two kinds of reactions are present in ReaxFF: valence
interactions and non-covalent (non-bonded) interactions. Valence interactions result from
overlap of electron wave functions. Non-covalent interactions consist of van der Waals
and Coulomb interactions, which account for repulsion at short distances from Pauli
exclusion principle and attraction at long distances due to dispersion.
In ReaxFF, van der Waals and Coulomb forces are screened by a taper function.
The taper function is a distance-dependent 7
th
order polynomial,
(2.61)
where Tap
7
= 20, Tap
6
= -70, Tap
5
= 84, Tap
4
= -35, Tap
3
= 0, Tap
2
= 0, Tap
1
= 0, Tap
0
=
1, and r
ij
is the distance between atoms i and j. Non-covalent energies are calculated for
all atom pairs separated by less than a cutoff length r
cnb
. When calculating a non-covalent
energy and force between a pair of atoms i and j, the energy and the force are multiplied
by the above taper function. The coefficients in the taper function are chosen such that all
1
st
, 2
nd
, and 3
rd
order derivates with respect to r
ij
are continuous and vanish at the cutoff
boundary.
33
In ReaxFF, van der Waals interaction is accounted for by a shielded Morse
potential given by
, (2.62)
(2.63)
Using a shielded potential prevents from very large repulsions between bonded atoms and
also between atoms sharing a valence angle. Coulomb interactions are also calculated by
a shielded potential as
(2.64)
€
H(r
ij
) = J
i
δ
ij
+
Tap(r
ij
)
r
ij
3
+(1γ
ij
)
3
[ ]
1 3
, (2.65)
€
H(r
ij
) = J
i
δ
ij
+
Tap(r
ij
)
r
ij
3
+(1γ
ij
)
3
[ ]
1 3
. (2.66)
All parameters in 2-body energy functions are listed in Table 2.6
34
Table 2.6: Parameters in 2-body energy functions
σ-, π-, and π π-bond energy parameters
bonding energy parameters
van der Waals energy parameter
van der Waals parameter
van der Waals parameter
van der Waals shielding
van der Waals parameter
Coulomb parameter
2.2.3.3 3-body Interactions
Valence angle, penalty, 3-body conjugation, and hydrogen bonding energies are
defined for triplets of atoms. Valence angle refers to an angle formed between three
atoms across two bonds, one atom being common between the two bonds; for example,
valence angle in water molecule is about 104.5˚. Equilibrium valence angle depends
on hybridization, over/under coordination, and number of lone pair electrons of central
atom in the triplet. The hybridization of the central atom, itself, depends on sum of π-
bond orders (SBO) of the central atom. In the absence of over/under coordination and
lone pair electrons, equilibrium valence angle changes from 109.5˚ for sp
3
hybridization
(0 π-bond) to 120˚ for sp
2
(1 π-bond) to 180˚ for sp (2 π-bonds). In ReaxFF, is defined
as
€
Θ
0
(SBO) =π −Θ
0,0
1−exp −p
val10
2−SBO2 ( )
[ ] { }
, (2.67)
(2.68)
35
€
SBO
j
= BO
jm
π
+BO
jm
ππ
( )
m∈n( j)
∑
+ 1− exp −BO
jm
8
( )
m∈n( j)
∏
⎡
⎣
⎢
⎢
⎤
⎦
⎥
⎥
−Δ
j
angle
− p
val8
n
lp,j
( )
, (2.69)
(2.70)
All parameters in equations of this section are listed in Table 2.7. Energy contribution
arising from deviation of valence angle Θ
ijk
from its equilibrium value is given by
(2.71)
€
f
7
(BO
ij
) =1−exp −p
val3
BO
ij
p
val4
( )
, (2.72)
€
f
8
(Δ
j
) = p
val5
− p
val5
−1
( )
2 +exp p
val6
Δ
j
angle
( )
1+exp p
val6
Δ
j
angle
( )
+exp −p
val7
Δ
j
angle
( )
. (2.73)
The function E
val
is defined such that energy contribution from a valence angle
approaches zero when bonds forming the valence angle dissociate.
When a valence angle is formed by two double bonds, such as propadiene,
deviations from the equilibrium angle is penalized with an additional energy term
(2.74)
(2.75)
In chemistry, conjugation refers to delocalization of π-bonds in a group, which
lowers the overall energy of the group and increases its stability. In ReaxFF, there are
two terms accounting for conjugation. 4-body conjugation term is devised for
36
incorporating the stability of conjugated hydrocarbon systems such as benzene while 3-
body conjugation term takes into account the stability of groups such as nitro group (—
NO
2
). 3-body conjugation energy E
coa
is calculated as
(2.76)
Normally, a hydrogen bond refers to attractive interaction between a hydrogen
atom of a molecule/group with a fluorine, oxygen, or nitrogen atom of another
molecule/group. However, in ReaxFF, a hydrogen bonding energy is associated with any
bonded triplet in which the central atom is hydrogen, and is defined as
€
E
hbond
= p
hb1
1−exp p
hb2
BO
ij
( )
[ ]
exp p
hb3
r
hb
0
r
jk
+
r
jk
r
hb
0
−2
⎛
⎝
⎜
⎞
⎠
⎟
⎡
⎣
⎢
⎢
⎤
⎦
⎥
⎥
sin
8
Θ
ijk
2
⎛
⎝
⎜
⎞
⎠
⎟
k
∑
j
∑
.
i
∑
(2.77)
37
Table 2.7: Parameters in 3-body energy functions
valence angle energy parameter
valence angle parameters
valency/lone pair parameter
, valence angle parameters
equilibrium valence angle parameter
atom valency
penalty energy parameter
double bond/angle parameter
, double bond/angle parameters for overcoordination
3-body conjugation energy parameter
valency angle conjugation parameters
hydrogen bond energy parameter
, hydrogen bond parameters
hydrogen bond radius
2.2.3.4 4-body Interactions
Torsion angle and 4-body conjugation energies are defined for quartet of atoms.
For each four atoms connected consecutively, the angle between the planes formed by the
first three and the last three atoms is called torsion or dihedral angle. More precisely, two
triplet of atoms (i, j, k) and (j, k, l) define a dihedral or torsion angle ω
ijkl
as
€
w
ijkl
= cos
−1
( ˆ n
ijk
⋅ ˆ n
jkl
), (2.78)
where and are unit vectors normal to planes defined by (i, j, k) and (j, k, l) triplets,
respectively. In ReaxFF, torsion angle energy contribution E
tors
is introduced to ensure
that dihedral angles remain close to their equilibrium values.
38
, (2.79)
€
f
10
(BO
ij
,BO
jk
,BO
kl
) =
1−exp −p
tor2
BO
ij ( ) [ ]
1−exp −p
tor2
BO
jk ( ) [ ]
1−exp −p
tor2
BO
kl
( ) [ ]
, (2.80)
€
f
11
(Δ
j
,Δ
k
) =
2 +exp −p
tor3
Δ
j
angle
+Δ
k
angle
( ) [ ]
1+exp −p
tor3
Δ
j
angle
+Δ
k
angle
( ) [ ]
+exp p
tor4
Δ
j
angle
+Δ
k
angle
( ) [ ]
. (2.81)
Torsion angle energy contribution vanishes when one of valence angles Θ
ijk
or Θ
jkl
approaches π. All parameters in equations of this section are listed in Table 2.8.
As mentioned above, 4-body conjugation energy accounts for stability of
hydrocarbon conjugated-systems. The 4-body conjugation energy E
conj
is defined as
€
E
conj
= p
cot1
f
12
(BO
ij
,BO
jk
,BO
kl
) 1+ cos
2
ω
ijkl
−1
( )
sinΘ
ijk
sinΘ
jkl [ ]
l
∑
k
∑
j
∑
i
∑
, (2.82)
€
f
12
(BO
ij
,BO
jk
,BO
kl
) =
exp−p
cot 2
BO
ij
−1.5
( )
2
[ ]
exp −p
cot 2
BO
jk
−1.5
( )
2
[ ]
exp −p
cot 2
BO
kl
−1.5
( )
2
[ ]
.
(2.83)
When consecutive atoms have bond orders of 1.5 (such as benzene or other aromatics),
contribution of 4-body conjugation energy is the largest.
39
Table 2.8: Parameters in 4-body energy functions
torsion energy parameters
torsion parameter
torsion/BO parameter
, torsion overcoordination parameters
4-body conjugation energy parameter
4-body conjugation parameter
2.3 Parallel Molecular Dynamics Simulations
Nowadays, scientific computing extensively relies on usage of computer clusters
consisting of thousands of CPU processors, GPU processors, and accelerators. There are
two reasons for which one may want to use more than one processor for a simulation:
simulating a larger system or simulating the system in hand faster. For a reason, which
will be clear shortly, the latter is more challenging. In both cases, the main task for
designing a parallel scheme is to divide computation load of the simulation between
processors. In MD simulations, there are different approaches for dividing the
computation load among which spatial decomposition and force decomposition are the
most common.
Dividing the computation load between processors is associated with two
overheads: communication and load balancing. The former is intrinsic to any parallel
scheme; there is always a communication overhead even for a homogeneous system. By
increasing number of processors, the communication overhead increases, and at some
point compensates the benefit of less computation load per processor. Therefore, for each
simulation, there exist an optimal number of processors, which produces the minimum
40
simulation time. Considering the communication overhead, it is clear that when the
number of particles per processor is small, the ratio between computation and
communication loads decreases resulting in a smaller optimal number of processors and
consequently more contingent limitation for making a simulation shorter.
In inhomogeneous systems, unbalanced load is another difficulty arising in
parallel computing. In general, there are two approaches for load balancing: static and
dynamic. In the former, load balancing takes place only once at the beginning of the
simulation while in the latter it repeats after a number of time steps. Dynamic load
balancing inevitably results in extra communication load, since after a number of time
steps the load of each processor should be evaluated, and some load should be migrated
between processors both of which require communication.
Recent processors have several cores. Multi-core processors are basically
different processors sharing their memory. Cores with a shared memory can
communicate faster as compared with processors having separate memories. There are
three common parallel programming libraries: Message Passing Interface (MPI),
OpenMP and Parallel Virtual Machine (PVM). MPI is a library for distributed-memory
parallel programming; OpenMP is a library for shared-memory parallel programming;
PVM is a library to use both shared- and distributed-memory architectures
simultaneously.
Although intra-processor communication is faster, but as of present, overhead of
initializing OpenMP threads is large such that using MPI alone usually gives a better
performance. For this reason, our parallel code used for simulations presented in this
41
dissertation solely use MPI library. Moreover, our systems are not extensively
inhomogeneous, so we do not utilize any load balancing scheme and load balancing will
not be discussed in the following.
2.3.1 Spatial Decomposition
As mentioned above, the main task in designing a parallel scheme is to divide the
computational load of the simulation between processors such that the simulation is as
fast as possible. To this end, communication and unbalanced load overheads should be
considered. There are numerous approaches for partitioning of the computational load
some of which are application-specific. Since designing and implementing a MD
program requires a considerable amount of effort and time, having a general-purpose MD
program is a wise approach—except for specific projects of a great importance—and
spatial decomposition is suitable for general-purpose MD programs intended to be used
for simulations with a large number of particles.
In the spatial decomposition approach, the simulation box is divided into domains
and each domain is assigned to one processor. Each processor is responsible for
computations needed for all atoms in its domain. We used MPI library for
communication between processors and domains in our program are of equal size. Let us
assume that P processors are available, and the MD box is partitioned into P = P
x
× P
y
×
P
z
domains, where P
x
, P
y
, and P
z
represent number of partitions along x, y, and z
directions, respectively. MPI assigns a serial ID (sid) to each processor. We assign a
42
vector ID (vid) to each processor, which is used to determine the processor to which an
atom belongs. Relation between sid and vid is established as
(2.84)
Fig. 2.5. depicts a specific case of the above relation with P
x
= 4, P
y
= 3, and P
z
= 2. The
coordinates of each atom determine the vid of the processor to which it belongs. Then,
the above relation provides the sid of that processor.
Fig. 2.5: Schematic illustration of correspondence between vid and sid. L
x
, L
y
, and L
z
denote lengths of the MD box in x, y, and z directions, respectively. These lengths will be
used to determine the vid of the processor to which an atom belongs and then the relation
between serial and vector IDs will determine the sid of the processor. Picture courtesy of
Aiichiro Nakano.
2.3.2 Communication Pattern
In the spatial decomposition scheme, each processor computes the trajectories of
the particles in the domain to which it is assigned. By being in charge of a particle, we
mean keeping all information of the particle such as coordinates and velocities, and also
performing all computation needed for that particle such as potential and force
43
calculations. Particles located close to boundaries of a domain interact with particles in
neighboring domains; therefore, in order to calculate potential and force, at each MD step,
each processor needs to obtain coordinates of atoms located in the neighboring domains
within a couple of cutoff distances (the distance depends on natures of interactions
present in the simulation e.g., 1 cutoff distance is sufficient if only pairwise interactions
are present).
Communication overhead depends on volume of the data and number of messages
to be transferred as well as pattern of the communication. The dependence of
communication overhead on volume of the data to be transferred is obvious because of
limited bandwidth of all communication apparatus. The dependence on number of
messages to be transferred arises from the latency associated with each message. This
factor can become significant if processors to be used are located far away from each
other say in different countries. Communication overhead also depends on the pattern of
the communication. In current supercomputers, it is common that each processor has
direct communication links to its up, down, right, left, front and back neighbors in the
cluster. If so, communication of each processor with these 6 neighbors is considerably
faster than with other processors in the cluster.
Considering the factors affecting communication load and also considering the
goal of having a general-purpose program, we used a communication pattern in which
each processor communicates with its up, down, left, right, front and back neighbor
processors. The communication operation consists of two functions: caching and
migration. To cache required coordinates of particles from 26 neighbor processors, 6
44
messages are transferred. First, each domain is extended in one direction say x. Then the
extended domain is extended in another direction say y, and finally the new extended
domain is extended in the remaining direction say z as shown in Fig. 2.6 for a two
dimensional case.
Fig. 2.6: Schematic of communication pattern in two dimensions. First, each domain is
extended from left and right. Then, the extended domain is extended again from up and
down directions. Picture courtesy of Aiichiro Nakano.
After each MD step, some particles leave their domain. The migration function
passes the information of these particles to their new processors via 6 messages the same
way as explained for cache function. Using this communication pattern, communicating
with 26 neighbors is achieved by 6 messages resulting in a low latency load. Moreover,
this pattern is suitable for architectures having direct links between a processor and its up,
down, left, right, front and back neighbor processors, which is quite common in
supercomputers as mentioned. We note that the spatial decomposition scheme and the
communication pattern discussed here, implicitly, assume that interactions present in the
simulation have a cutoff distance. In the case of long-range interactions, other methods
such as fast multiple method (Greengard and Rokhlin 1997) and Ewald summation
45
(Frenkel and Smit 2002) can be used to reduce the complexity of force calculation from
O(N
2
) to O(N) (N being the number of particles in the simulation box). Accordingly,
appropriate communication patterns should be adopted as well.
2.3.3 Performance of Parallel Programs
Let us consider a problem of size W
P
(computation load equal to W
P
) to be solved
on P parallel processors using a specific parallel program. Let be the execution
time required to accomplish this task. Speed of the program is defined as
. Speedup using P processors is defined as the ratio between the
speed of P processors and that of 1 processor i.e., . The
ideal speedup on P processors is P; therefore, parallel efficiency using P processors is
defined as .
Now, let us consider a problem of size W. As mentioned earlier, there are two
reasons for using parallel computers: solving a larger problem or solving the same
problem faster. Accordingly, there are two scalability measures. In weak (isogranular)
scalability test, the work per processor is constant. Therefore,
(2.85)
(2.86)
In a strong scalability test, the total work is constant. Hence,
46
(2.87)
(2.88)
Let us consider a problem in which a fraction f (0 < f < 1) of the work is inherently
sequential and cannot be parallelized. In this case,
(2.89)
This result known as Amdahl’s law reveals that the upper bound for speedup (for strong
scaling) is determined by f. For example, let us consider a problem with
€
f =0.01. The
maximum speedup for this problem is
€
1/ f =100; therefore, using more than O(100)
processors for such a problem is not efficient.
47
3 Structure and Dynamics of Shock-induced
Nanobubble Collapse in Water
3.1 Simulation Setup
We have performed MD simulations of shock-induced collapse of nanobubbles in
water. The MD approach is well suited to study this problem, because it can provide
direct information about the structure and dynamics of nanobubble collapse over
microscopic spatio-temporal scales. Contained in this microscopic information are some
of the subtle but salient features of bubble collapse that experiments or continuum
simulations may not be able to capture.
The MD simulations reported here are based on quantum-mechanically informed
reactive force fields ReaxFF, which can accurately describe bond dissociation/formation
and chemical reactions in the system (van Duin, Dasgupta et al. 2001). In ReaxFF, charge
transfer between atoms changes dynamically with changes in the atomic configuration
and atomic charges are obtained by minimizing the electrostatic energy at every MD step
(Rappe and Goddard 1991). The ReaxFF also includes three- and four-body interactions.
The ReaxFF potential for water was trained against a wide range of quantum mechanical
data, including binding energies for n-water clusters (n = 2-30), proton transfer reactions
in H
3
O
+
-H
2
O and OH-H
2
O dimers, concerted proton transfer reactions in [H
2
O]
n=2-6
clusters and dissociation energies for H-H, O-H, O-O and O=O bonds. The parameters
were validated for bulk water by comparing the density, cohesive energy, self-diffusion
48
coefficients and O-H, O-O and H-H radial distribution functions with experiments.
Using scalable fast reactive force field algorithm F-ReaxFF (Nakano, Kalia et al.
2007; Nomura, Kalia et al. 2008), we have performed MD simulations of water subjected
to planar shock with and without a nanobubble. Initial dimensions of the MD box are
38.5 nm, 18.6 nm, and 16.6 nm in the x, y, and z directions, respectively, and the system
contains a million atoms. We equilibrate the system in the (N, V, T) ensemble i.e.,
constant number of atoms, constant volume, and constant temperature with PBC in all
directions. The initial mass density and temperature are 0.98 g/cc and 300 K, respectively.
After equilibration, we insert nanobubbles of different diameters (D = 6, 8, and 10 nm) at
the center of the MD cell by removing 90% of water molecules from spheres of
corresponding diameters; see Fig. 3.1. Although the nanobubbles are thermodynamically
unstable, they are subjected to shock compression so rapidly (in about a pico second) that
they do not have time to shrink on their own over such a short duration.
49
Fig. 3.1: Schematic of the simulation cell. The grey plate is the momentum mirror.
These equilibrated systems are then subjected to planar shocks with particle
velocities u
p
= 0.4, 1.0, 1.5, 2.0, 2.5, 3.0, or 3.5 km/s. Shocks are applied in the x
direction using a momentum mirror positioned at the left end of the box. At time t = 0 all
atoms are given an initial particle velocity, u
p
, towards the momentum mirror, and
whenever an atom crosses the y-z plane of the mirror, the x-component of its velocity is
reversed. In shock simulations, we insert a 2-nm-thick vacuum layer at the end of the MD
cell in the x direction (shock direction), turn off the thermostat coupling, and apply PBC
in lateral directions to minimize surface effects normal to the shock direction. The
distance between a nanobubble and its nearest mirror image (≥ 6.6 nm) is much larger
than the cut-off length (1 nm) of ReaxFF, which rules out any interaction between the
nanobubble and its images. The integration time step in these MD simulations is 0.1 fs.
50
3.2 Structure of Water under Shock
We first performed shock simulations without the nanobubble to validate the
reactive force field of water. Fig. 3.2 shows the setup for this simulation in the inset, and
MD (red circles) and experimental (blue crosses) (Rybakov and Rybakov 1995) results
for the shock velocity, u
s
, as a function of the particle velocity, u
p
, i.e. the Hugoniot curve.
Shock velocity is obtained from the difference in the shock-front boundaries at two time
frames. In each frame, the abrupt change in the density identifies the location of the
shock front. The simulation results for the Hugoniot curve are in good agreement with the
experimental results.
Shock produces significant structural changes in water. Fig. 3.3 shows that the
oxygen-oxygen radial distribution function (RDF) for water under planar shock is
significantly different from that for ambient water. The maxima and minima in the RDF
for shocked water (u
p
= 3 km/s) are sharper and shifted with respect to those in ambient
water. The most interesting structural changes occur at u
p
= 1 km/s. In the primary shock
region, we find molecular clusters in which a water molecule has 8 nearest neighbors in a
body-centered cubic lattice configuration, indicating the nucleation of ice VII; see Fig.
3.4. In this ice VII structure, the central oxygen atom is connected to four of its nearest-
neighbor oxygen atoms through hydrogen bonds (Jorgensen and Worlton 1985; Klotz,
Bove et al. 2009). We do not find an ice-like structure at particle velocities u
p
= 2.5, 3.0
and 3.5 km/s, which is in agreement with experimental observations of ice VII only for
particle velocities between 0.75 and 2 km/s (Rybakov and Rybakov 1995)
51
Fig. 3.2: Experimental (blue crosses) and MD (red circles) Hugoniot compression curves
for particle velocities between 0.4 and 3 km/s. The inset shows the simulation cell and the
momentum mirror (grey plate).
52
Fig. 3.3: Oxygen-oxygen pair-correlation functions, g(r), for ambient and shocked water
with and without a nanobubble as a function of distance at u
p
= 3 km/s. The nanobubble
diameter is 10 nm.
53
Fig. 3.4: Snapshot of a water molecule cluster that forms an ice-VII-like structure in the
compressed region at u
p
= 1 km/s. Red and white spheres represent oxygen and hydrogen
atoms, respectively. Dotted lines indicate the directions of nearest neighbor oxygen atoms
(within 3Å from the central oxygen atom).
3.3 Dynamics of Shock-induced Nanobubble Collapse
Fig. 3.5(a) shows how the nanobubble volume (normalized by the initial
nanobubble volume) and shape of the largest bubble change with time as a result of shock
compression. The insets show oxygen atoms of water molecules at the periphery of the
largest nanobubble (D = 10 nm) before the shock wave strikes the nanobubble (t = 1.7 ps)
and just before the nanobubble collapses (t = 3.3 ps). The volume versus time data
correspond to nanobubble diameters D = 6, 8 and 10 nm and particle velocity u
p
= 3 km/s.
The bubble collapse time τ is estimated to be 1.1, 1.4 and 1.7 ps for D = 6, 8 and 10 nm,
54
respectively. According to the Rayleigh formula ( , where ρ is the mass
density and ΔP is the pressure difference across the bubble surface), τ for the three
nanobubble sizes are 0.8, 1.1 and 1.4 ps. The differences between the MD results and
Rayleigh formula arise for the following reasons: (1) in Rayleigh collapse it is assumed
that the bubble is within a uniform fluid, whereas the pressure and density around the
nanobubble become non-uniform when the shock front reaches the nanobubble; and (2)
unlike MD simulations, Rayleigh equation does not include viscosity and surface tension
effects arising from interatomic interactions.
It should be pointed out that not all nanobubbles collapse under shock.
Experimentally, surface gas nanobubbles are stable under tensile and compression
pressures of 6 MPa (Borkent, Dammer et al. 2007). Here, we are dealing with cavitation
nanobubbles subjected to compression waves with high pressures ranging between 2-19
GPa (depending on the particle velocity), which is why nanobubbles collapse (see Fig.
3.6 for more details).
Fig. 3.5(b) and 3.5(c) show velocity profiles of water molecules around a
shrinking nanobubble (D = 10 nm) at u
p
= 3 km/s. At t = 3.3 ps the nanobubble has
shrunk significantly but not collapsed completely under shock compression. The snapshot
in Fig. 3.5(c) is taken immediately after the nanobubble collapses (t = 3.9 ps). Here,
contours of the magnitude of velocity are color-coded and each arrow represents the
direction and the arrow color the magnitude of the average molecular velocity (averaged
over all molecules in a voxel of length 0.5 nm) at that position. Regions inside the white
55
dash lines are enlarged in the insets, and the white arrows indicate the direction of shock
propagation. We observe that water molecules around the top and bottom of the
nanobubble change directions and their average molecular velocities point towards the
center of the nanobubble soon after the shock wave reaches the proximal side of the
nanobubble. This focusing feature of high-speed molecules from the onset of shrinkage to
complete collapse of a nanobubble is akin to the jetting phenomenon observed
experimentally in microbubble collapse. The difference between our simulation and
experiment (Ohl and Ikink 2003) is that nanojets do not elongate as much as microjets
beyond the distal side of collapsed bubbles. The formation of nanojets has also been
reported in an MD simulation of pressurized fluid injection through a nanoscale nozzle
(Moseler and Landman 2000). Fig. 3.5(c) shows molecular motion around the collapsed
nanobubble and lateral flow at the periphery of the collapsed region. The lateral flow is
smaller than that observed in the vicinity of a solid surface or a lipid bilayer (Ohl, Arora
et al. 2006; Ohl, Arora et al. 2006).
56
Fig. 3.5: (a) Normalized nanobubble volume, Ω(t), versus time, t, for nanobubbles of
initial sizes D = 6, 8 and 10 nm. Insets show snapshots of the largest nanobubble before
the shock wave strikes (t = 1.7 ps) and before the nanobubble collapses (t = 3.3 ps). (b)
and (c) show velocity vector fields in a cross section of the system just before (t = 3.3 ps)
and soon after (t = 3.9 ps) the largest nanobubble collapses. The velocities are measured
relative to the reference frame of momentum mirror. In (b) the region inside the white
dash lines is enlarged to show a focused nanojet in the inset; in (c) the lateral flow
indicates the onset of nanojet disintegration.
57
Fig. 3.6: Time evolution of the shape of the nanobubble of diameter D = 10 nm at u
p
= 3
km/s.
At the onset of nanobubble collapse, we observe a sudden increase in the
translational kinetic energy of molecules at the shock front and also an appreciable
increase in the rotational energy of those molecules. In the final stages of nanobubble
collapse, the molecular vibrational energy also increases to the same level as the
rotational energy. For u
p
= 3 km/s and D = 10 nm, the maximum rotational and
vibrational energies are about 0.2 eV, which is well below the self-dissociation energy
[2H
2
O → OH
-
+ H
3
O
+
] in pure water; see Fig. 3.7. However, we observe significant
number of positive and negative ions produced by nanobubble collapse. This issue will be
discussed in more detail in Chapters 4 and 6.
58
Fig. 3.7: Rotational (a) and vibrational (b) temperature profile along the direction of
shock propagation—x. Rise in vibrational temperature lacks the one in rotational
temperature. The rise in rotational temperature occurs at t = 2.9 ps while the rise in
vibrational temperature takes place at t = 3.4 ps.
Fig. 3.8 shows effects of the shock amplitude (i.e., particle velocity) and
nanobubble size on velocity vector fields of water molecules in and around the
nanobubble just after it collapses. Panels (a), (b) and (c) display velocity vector fields for
a fixed particle velocity, u
p
= 3 km/s, but different initial diameters D of the nanobubble.
Evidently, the size of the nanojet (ellipsoidal red region) increases with an increase in the
nanobubble diameter. Panels (d), (e) and (f) show the effect of the particle velocity on the
molecular velocity vector fields in systems containing the same size nanobubble (D = 10
nm). The high-speed region does not grow as much with an increase in the particle
velocity as with the nanobubble size.
59
Fig. 3.8: Effects of particle velocity and nanobubble size on nanobubble collapse.
Snapshots in panels (a)-(c) are taken for u
p
= 3 km/s and nanobubble diameter D = 6, 8
and 10 nm, respectively. Snapshots in panels (d)-(f) are for a fixed nanobubble size (10
nm) and u
p
= 2.5, 3.0 and 3.5 km/s, respectively. The size of the high-speed region
(ellipsoidal red regions) increases more with an increase in D than with u
p
. All snapshots
are taken immediately after the nanobubbles collapse.
From the onset of nanojet formation and disintegration, we have determined the
penetration length l
jet
and persistence time τ
jet
for the jet to examine how these quantities
are affected by the particle velocity u
p
and initial nanobubble radius R. For all particle
velocities we find l
jet
increases linearly with R. For example, we find that l
jet
increases
from 7.8 nm to 14.7 nm as R increases from 3 nm to 5 nm. Experimentally, lengths of jet
tip in micron-to-millimeter-sized bubbles also scales linearly with R (Kodama and
Tomita 2000; Ohl and Ikink 2003). This agreement is somewhat surprising, because
viscosity and surface tension may play larger roles in nanobubbles than in micron-to-
millimeter-sized bubbles. We find the nanojet persistence time is longer than the
nanobubble collapse time by about 0.2 ps. Also the nanojet velocity, V
jet
= l
jet
/τ
jet
, is
60
higher than the shock velocity u
s
(and hence also u
p
), and V
jet
increases with an increase
in the nanobubble size.
Fig. 3.9 shows the local pressure from the virial expression (Sun, Wang et al.
2006) for a nanobubble of diameter D = 10 nm at u
p
= 3 km/s. The pressure in the
compressed region remains constant and the wave front remains planar until the shock
wave reaches the nanobubble; see Fig. 3.9(a). The pressure (18.8 GPa) is in reasonable
accord with the estimate (19.5 GPa) from the jump condition, (Warsi
1998). Fig. 3.9(b) shows the wave front becomes concave due to pressure gradient in the
compressed region during the nanobubble collapse. As the nanobubble shrinks, the
pressure at the shock front decreases. The shock front becomes planar again just before
the complete collapse of the nanobubble; see Fig. 3.9(c). When the primary shock wave
hits the proximal side of the nanobubble, surrounding water molecules move into the
nanobubble with high speeds. On reaching the distal side of the nanobubble, these high-
speed molecules give rise to a secondary water hammer shock wave with a maximum
pressure of 29 GPa [Fig. 3.9(d)] while the primary shock continues to move forward with
a pressure of 15 GPa. The water hammer shock propagates backward (opposite to the
primary shock), spreading spherically with a velocity of 8 km/s as the pressure decreases;
see Fig. 3.9(e) and (f). Snapshots in Fig. 5 show that the water hammer shock is not close
to the lateral boundaries and so the interaction between mirror-image waves is negligible.
Despite the differences in length scale and pulse duration, significant pressure
amplification due to water hammer shock wave and its rapid dissipation have also been
61
observed in experiments (Brujan, Keen et al. 2002) and continuum simulations (Johnsen
and Colonius 2008).
Fig. 3.9: [(a)-(c)] show the primary shock wave at t = 1.6 ps, 2.9 ps, and 3.3 ps,
respectively. Here, the particle velocity, u
p
= 3 km/s and D = 10 nm. [(d)-(e)] show the
water hammer shock at t = 3.5 ps, 3.7 ps, and 3.9 ps, respectively. The color code
corresponds to pressure and shock waves are identified by pressure discontinuities. The
length scale bar for these panels is shown in (a).
62
4 Poration of Lipid Bilayers by Shock-induced
Nanobubble Collapse
4.1 Simulation Setup
We have performed MD simulations to study the impact of shock waves on
nanobubbles in the vicinity of a dipalmitoylphosphatidylcholine (DPPC) phospholipid
bilayer embedded in water. We use the simple point charge (SPC) model for water
(Berendsen, Postma et al. 1981; Berweger, Vangunsteren et al. 1995). Our simulations
reveal molecular mechanisms of nanobubble shrinkage and collapse by shock waves, the
nanojet formation and penetration into lipid bilayers, and the resulting changes in
structural and dynamical properties of bilayers including their recovery after the passage
of shock waves.
First, we checked the validity of the SPC model for water under shock. It is
reasonable to use the simple point charge (SPC) model for water under shock even
though the model does not allow bond breaking. We have checked this by performing
MD simulations for shock-induced nanobubble collapse in water using a reactive force
field (ReaxFF), which can accurately describe bond breaking/formation and chemical
reactions in the system. An MD simulation was performed to equilibrate the water system
by itself for 1 ns using a time step of 2 fs. PBC were imposed during equilibration. To
apply shock, we inserted a vacuum layer of thickness 2 nm at the end of the MD box in
the x direction and moved the system with a constant particle velocity u
p
towards a
63
momentum mirror (Nomura, Kalia et al. 2007) i.e., along negative x direction in the inset
of Fig. 4.1. The mirror reverses the x component of atomic velocity if an atom crosses its
plane. This generates a planar shock in the positive x direction. The shock velocity u
s
is
determined by monitoring the shock front i.e., discontinuity in pressure or mass density
of water at two instants of time. Fig. 4.1 shows that the MD results for u
s
as a function of
u
p
are in good agreement with experimental data (Rybakov and Rybakov 1995).
Next, we performed MD simulations to equilibrate initial configurations of lipid
bilayers and water molecules at temperature T
i
= 300 and 323 K and pressure P = 1 bar
using a time step of 2 fs. We use GROMACS (Hess, Kutzner et al. 2008) to simulate the
SPC model for water and DPPC bilayer. To create systems large enough to accommodate
nanobubbles of different diameters (D) (to be inserted later), we repeat the configuration
taken from (Tieleman and Berendsen 1996) in the y-z plane (where x is along the bilayer
normal and shock direction) 3×3, 5×5, 10×10 times for D = 10, 20, 40 nm, respectively.
The parameters for bonded and non-bonded interactions as well as partial charges on
atoms in the DPPC bilayer system are taken from (Tieleman and Berendsen 1996; Berger,
Edholm et al. 1997). The systems with bubbles of diameters D = 10, 20 and 40 nm
contain about 1.9, 6.2, and 30.8 million atoms, respectively, and the dimensions of the
corresponding MD cells are 19×19×61, 34×34×82, and 64×64×92 nm
3
.
64
Fig. 4.1: Shock velocity versus particle velocity. The simulation results for SPC water are
in good agreement with experimental data. The inset shows the setup for shock
simulation. The gray plate is the momentum mirror.
During equilibration the DPPC bilayer and water molecules are coupled to a
thermostat and barostat (Berendsen, Postma et al. 1984), bond lengths are kept constant
(Hess, Bekker et al. 1997), and PBC are imposed in all directions. After equilibration,
nanobubbles are inserted in the vicinity of the lipid bilayer by removing water molecules
within a sphere of diameter 10, 20 or 40 nm; see Fig. 4.2. The distance between the right
edge of the nanobubble and the left boundary of the lipid bilayer is 1.9 nm. Shock is
applied using a momentum mirror located at x = 0.2 nm. A vacuum layer of width 2 nm
is inserted on the right side of the MD cell before applying the shock. Pressure and
65
temperature couplings as well as bond-length constraints are turned off in shock
simulations. A number of simulations were performed for different nanobubble diameters
(D = 10, 20 and 40 nm), particle velocities (u
p
= 0.4, 0.7 and 1 km/s), and initial
temperatures (T
i
= 300 and 323 K) (Albon and Sturtevant 1978).
Fig. 4.2: Schematic of the MD box containing a nanobubble in water and a lipid bilayer.
The gray plate is the momentum mirror.
4.2 Effects of Nanobubble Collapse on Lipid Bilayers
When a planar shock front hits the proximal side of a nanobubble, water molecules
from the nanobubble periphery accelerate towards the center of the nanobubble and form
a nanojet. The size of the nanojet depends on the particle velocity and nanobubble
diameter. As the particle velocity increases from 0.4 km/s to 1.0 km/s, the number of
66
water molecules in the nanojet for a fixed value of D increases by an order of magnitude.
Fig. 4.3 displays instantaneous molecular velocities averaged in voxels of dimension 0.5
nm for a nanobubble of initial diameter D = 40 nm under the impact of a shock front
moving with velocity u
p
= 0.7 km/s. Fig. 4.3(a) shows that velocities of water molecules
in the domain of the shrinking nanobubble are focused in the form of a nanojet. As the
particle velocity increases from 0.4 km/s to 1.0 km/s, the average x-component of
molecular velocities inside the nanojet increases from 2.6 km/s to 3.5 km/s for all
simulated nanobubble sizes. For D = 40 nm, we find the length of the nanojet l
jet
is 57 nm.
Our results for other nanobubble sizes and particle velocities also indicate that the length
of the jet scales as l
jet
≈ 1.5D. Surprisingly, the same linear scaling has been observed in
experimental studies of shock-induced collapse of micron-to-millimeter-sized bubbles
(Kodama and Tomita 2000; Ohl and Ikink 2003). In Fig. 4.3(b) we show the interaction
between water molecules in the nanojet and the DPPC molecules in the bilayer. We
observe vortices in the collapsed bubble when water molecules bouncing back from the
bilayer encounter other water molecules in the incoming shock wave.
67
Fig. 4.3: Snapshots of velocity profile for the system with D = 40 nm, T
i
= 300 K and u
p
=
0.7 km/s. Arrows show the direction of average molecular velocities and the velocity
magnitudes are color-coded. Panel (a) shows a nanojet in the system at t = 20 ps. The
white vertical region is the bilayer. Panel (b) shows a spreading flow at t = 24 ps resulting
from the impact of the nanojet on the lipid bilayer.
Water molecules in the nanojet form a spreading flow after hitting the first leaflet of
the DPPC bilayer. Fig. 4.4 shows the averaged lateral velocity of water molecules in the
vicinity of a lipid bilayer versus the distance from the center of the bilayer for a
nanobubble of diameter 40 nm. The inset shows the results for a nanobubble of diameter
10 nm. In both systems u
p
= 0.7 km/s and T
i
= 300 K. The peak in the lateral-flow
velocity appears when the nanojet hits the bilayer (t = 21 ps for the larger nanobubble and
t = 7 ps for the smaller one). The distance over which the lateral velocity is larger than
the thermal velocity is half the nanobubble radius. Experiments on millimeter-sized
bubbles in the vicinity of a hard surface indicate that the distance is of the order of the
bubble radius (Ohl, Arora et al. 2006)). The differences between experimental and our
68
MD results are due to the fact that bubble sizes differ by several orders of magnitude and
the surfaces are soft in MD simulation and hard in experiments.
Fig. 4.4: Average lateral velocity of water molecules in the vicinity (1nm) of a lipid
bilayer versus the distance from the center of the bilayer for a nanobubble of diameter 40
nm. The inset shows the results for a nanobubble of diameter 10 nm. In both systems u
p
=
0.7 km/s and T
i
=300 K.
Fig. 4.5(a) shows the water density around the DPPC bilayer (blue region) just
before the bilayer is hit by the water nanojet from a collapsing nanobubble of diameter 40
69
nm at u
p
= 0.7 km/s. The curved blue region to the left of the bilayer indicates that the
nanobubble has not collapsed completely and the water density around the nanobubble is
close to the normal density of water. After the nanojet impact, the DPPC bilayer deforms
and becomes significantly disordered (see below). The water density around the bilayer
leaflet closer to the collapsed bubble increases to 1.5 g/cc. Fig. 4.5(b) shows that the
deformed bilayer is hemi-spherical. In addition, we observe water hammer shock when
water molecules in the nanojet hit the distal side of the nanobubble. This secondary water
hammer shock spreads spherically and its initial speed (first 4 ps after formation) is
approximately 1.6 km/s. The amplitude of the secondary shock decreases, but its velocity
increases with time.
Fig. 4.5: (a) and (b) are snapshots of the density of water at t = 20 and 28 ps. Here D = 40
nm, T
i
= 300 K and u
p
= 0.7 km/s. The central blue region is the lipid bilayer. Panel (a)
shows the nanojet traveling towards the distal side of the nanobubble. Panel (b) shows the
deformed bilayer and water hammer shock.
70
To characterize the evolution of disorder in the bilayer, we calculate the order
parameter tensor defined as (Tieleman and Berendsen 1996)
(5.1)
where is the angle between the ith “molecular vector” and the bilayer normal. The
“molecular vector” for the ith CH
2
unit has z component from C
i-1
to C
i+1
and the y
component in the plane containing C
i-1
, C
i
, and C
i+1
. The brackets denote an ensemble
average. From the diagonal elements of , the normalized order parameter at time t is
defined as
(5.2)
where t = 0 is the time at which the nanobubble is inserted in the system. Fig. 4.6 shows
how rapidly the impact of the nanojet partially destroys order in the lipid bilayer. After
impact, the order parameter within a radius of 5 nm from the center of the hemi-
spherically deformed bilayer decreases by 60% in about 4 ps. The disorder in the bilayer
is less farther from the deformation center. For example, the drop in the order parameter
20 nm from the center of the deformed bilayer region is 36%.
71
Fig. 4.6: Temporal changes in the normalized order parameter for the lipid bilayer region
directly impacted by the nanojet. At t = 16 ps, the density of water molecules around the
bilayer starts to increase, and hence the bilayer begins to show increasing disorder.
4.3 Poration and Recovery
The impact of the nanojet causes poration in the lipid bilayer. Fig. 4.7 shows
poration resulting from the impact of the collapsed nanobubble of initial diameter 40 nm
at u
p
= 0.7 km/s. The poration was calculated by dividing the impacted region of the
bilayer into pixels of size 0.1 nm and determining the area of empty pixels i.e., those
containing no lipid molecules. For the bilayer initially in the gel phase (Albon and
Sturtevant 1978) at room temperature (T
i
= 300 K), the nanojet impact increases poration
by a factor of 30 over its normal value before the nanojet impact; see Fig. 4.7(a). For the
bilayer initially in the liquid phase (Albon and Sturtevant 1978) at T
i
= 323 K the
72
poration increases by another factor of 5 relative to the poration in the gel phase; see Fig.
4.7(b). In the liquid phase at 323K, the maximum nanopore size is 0.7 nm as compared to
0.4 nm in the gel phase. The poration also depends on the particles velocity and
nanobubble diameter. At u
p
= 0.4 km/s, we do not observe any significant change in the
porosity of the gel phase for the three bubble sizes we have considered. However, at u
p
=
1 km/s, the maximum nanopore size increases to 0.3 nm for D = 10 nm, and it increases
linearly with the initial diameter of the nanobubble.
Fig. 4.7: Poration of lipid bilayers by collapsed nanobubbles. Here u
p
= 0.7 km/s and D =
40 nm. The bilayer was initially in the gel phase at T
i
= 300 K (a) and in the liquid phase
at T
i
= 323 K (b).
In the deformed DPPC bilayer that was initially in the liquid phase, the pores are
large enough (≥ 0.5 nm) to allow translocation of water molecules. Translocation events
are seen for u
p
= 1.0 km/s and D ≥ 10 nm and also for u
p
= 0.7 km/s and D = 40 nm.
Snapshots in Fig. 4.8 display translocation of a water molecule through a lipid bilayer
73
deformed by the nanojet impact from the collapse of a nanobubble of diameter D = 10 nm
at u
p
= 1.0 km/s. Here, translocation occurs because the nanojet impact creates a
considerable pressure difference (~ 9 GPa) across the bilayer. In Fig. 4.8, the first
snapshot on the left shows the bilayer and a water molecule (solid white circle) at time t =
4 ps, which is before the nanojet hits the bilayer. At t = 8 ps the nanojet impact deforms
and thins the bilayer in the middle, initiating the translocation of the water molecule. At t
= 12 ps the water molecule reaches the predominantly hydrophobic region in the middle
of the bilayer. The snapshot at t = 16 ps clearly shows the water molecule has moved
across the bilayer to the hydrophilic region of the second leaflet.
74
Fig. 4.8: Passage of a water molecule (white dot) across the lipid bilayer for the system
with D = 10 nm, T
i
= 323 K, and u
p
= 1.0 km/s. The water molecule, initially in the red
hydrophilic region on the left side of the bilayer, crosses the hydrophobic region (light
green) and ends up on the right side (red) of the bilayer.
We have also examined the healing of lipid bilayers and nanopores after the
passage of shock waves. Fig. 4.9 shows the relaxed configuration of the lipid bilayer in
which u
p
= 1.0 km/s, D = 10 nm, and T
i
= 300 K. Here, the MD simulation was performed
on a configuration of a deformed bilayer (including surrounding water molecules in an 11
nm wide region) embedded in water under ambient conditions. The simulation was run
for 5 ns in the microcanonical ensemble using PBC. After relaxation, the deformation in
the bilayer disappears.
75
Fig. 4.9: Deformed (a) and relaxed (b) lipid bilayer configurations. After relaxation, the
deformation in the bilayer disappears.
76
5 Mechano-chemical Pathways to H
2
O and CO
2
Splitting
5.1 Simulation Setup
We performed MD simulations to investigate chemical reactions triggered by
shock-induced collapse of a CO
2
-filled nanobubble in water. One interesting aspect of
this simulation is the sole existence of H
2
O and CO
2
molecules. The simulation is based
on a quantum-mechanically informed reactive force field (ReaxFF) to describe bond
formation and dissociation (van Duin, Dasgupta et al. 2001). We use a scalable parallel
algorithm (F-ReaxFF) to perform large reactive MD simulations on parallel computers
(Nakano, Kalia et al. 2007; Nomura, Kalia et al. 2008). We simulate shock-induced
nanobubble collapse in water without and with CO
2
molecules inside the nanobubble,
where the system without CO
2
molecules serves as a reference system. We start with a
MD cell filled with water molecules at a density of 0.98 g/cc and temperature T = 300 K.
The size of the MD cell containing N = 10
6
atoms is V = 35.8 × 18.6 × 16.6 nm
3
in the x,
y and z directions. After thermalization in the canonical ensemble (i.e., constant N, V, and
T) for 0.1 ps, the system is subjected to a planar shock with particle velocity (u
p
) 3 km/s
propagating in the x direction. At this point, the thermostat is turned off, a 2-nm-thick
vacuum layer is added to the end of the system in the x direction, and PBC along the
shock direction is removed. To apply the shock, all atoms are given a velocity of 3 km/s
in the negative x direction towards a momentum mirror located at x = 0. When an atom
77
passes across the momentum mirror, the x component of its velocity is reversed; see Fig.
5.1.
Fig. 5.1: Schematic of the primary system. All atoms are given a velocity u
p
in the
negative x direction towards a momentum mirror, and the x component of their velocities
is reversed when they pass the momentum mirror. The nanobubble with diameter 10 nm
contains 532 CO
2
molecules.
After 1.6 ps, the reference system is created by removing water molecules from a
sphere of diameter 10 nm at the center of the MD cell. At this moment (t = 0), the
distance between the shock front and the proximal point of the nanobubble is 2 nm. The
primary system is the same as the reference system except for the nanobubble being filled
with 532 CO
2
molecules. The CO
2
molecules are extracted from a CO
2
gas at T = 300 K
and positioned randomly inside the nanobubble. The MD time step in these simulations is
0.1 fs, and t = 0 corresponds to the moment at which the nanobubble is created.
78
5.2 Ionization of Water
The dynamics of shock-induced collapse of nanobubbles have been discussed in
Chapter 3. For the present simulations with u
p
= 3 km, the shock front velocity (u
s
) is 6.6
km/s. As the shock propagates towards the nanobubble, the pressure increases to 18.8
GPa in the compressed region behind the shock front. When the shock front strikes the
nanobubble, the surrounding water molecules move into the nanobubble with high speeds,
forming a nanojet. Molecules forming the nanojet have velocities greater than the shock
front and are focused towards the distal point of the nanobubble. Upon impact of the
nanojet on the distal side of the nanobubble, a secondary water hammer shock is
generated. While the primary shock continues to move forward, the water hammer shock
propagates spherically backward with a velocity (u
s
= 8 km/s) higher than the primary
shock velocity. When the water hammer shock forms, the local pressure rises to 29 GPa,
which is significantly greater than the pressure generated by the primary shock. The
pressure amplification dissipates rapidly as the water hammer shock spreads.
Fig. 5.2(a) shows the numbers of positive and negative ions as a function of time
in the reference system. Hydronium (H
3
O
+
) and hydroxyl (OH
-
) ions are the most
abundant positive and negative ions, respectively, though we also observe other positive
and negative ions such as H
+
or H
3
O
2
-
. At time t = 1.7 ps, when water hammer shock
forms, the number of ions starts to increase and reaches 190 at t = 2.5 ps. In the reference
system, there are equal numbers of positive and negative ions.
79
Fig. 5.2: The number of ions as a function of time in the reference system (a) and in the
primary system (b).
Short time and length scales involved in the autoionization of water preclude
direct experimental observations, but an ab initio MD simulation (Geissler, Dellago et al.
2001) revealed that rare fluctuations in solvation energies destabilizing O-H bonds,
followed by proton hopping through hydrogen bond “wires”, result in separated ions. If a
hydrogen bond wire connecting two ions separated by three or more neighbors breaks,
two stable ions may form. When the nanobubble collapses, high local pressure
accompanying the formation of the water hammer shock brings water molecules closer to
each other. The close distance between water molecules assists the destabilization of O-H
bonds and also strengthens hydrogen-bond networks, which together increase the
likelihood of autoionization events. As the water hammer shock spreads and local
pressure decreases, the number of ions stops increasing and fluctuates around 190.
Fig. 5.2(b) shows the numbers of positive and negative ions as a function of time
in the primary system and the spatial distribution of these ions are shown in Fig. 5.3. A
comparison between the reference and primary systems reveals that while numbers of
80
positive ions in both systems are approximately equal, there is an excess of negative ions
in the primary system. Moreover, the generation of ions starts earlier (at t = 1.4 ps). These
differences arise from H
2
O-CO
2
reactions in the primary system and will be discussed
further in the following. Ions and carbon atoms in the primary system at the end of the
simulation (t = 3.5 ps) are displayed in Fig. 5.4.
81
Fig. 5.3: Snapshots of the primary system showing spatial distribution of ions at (a) t =
0.1 ps before the shock front strikes the nanobubble, (b) t =1.7 ps when the water hammer
shock forms, and (c) t = 3.5 ps when the number of ions has reached a stationary value.
Positive and negative ions are shown in red and blue, respectively.
82
Fig. 5.4: Ions and carbon atoms in the primary system at t = 3.5 ps. (a) Side view. (b)
Front view. Positive and negative ions are shown in red and blue, respectively. Green
spheres depict carbon atoms.
5.3 Mechano-chemical Reactions between H
2
O and CO
2
Abundance of ions after the complete collapse of the nanobubble together with
proximity of H
2
O and CO
2
molecules under shock form an environment in which some
exotic pathways to H
2
O and CO
2
splitting become possible. Through close examination
of ionization events, we have identified the four dominant pathways summarized in Fig.
5.5. Each row in the figure represents one pathway, where the reactants and intermediate
states are shown respectively in the left and right panels:
1. In the first pathway in Fig. 5.5(a), a water molecule located between CO
2
and
other water molecules dissociates first. Then, the oxygen atom labeled 3 (hereafter
denoted as O(3)) and H(2) transfer to C(1), and the second hydrogen atom H(4) transfers
to O(5) of a neighboring water. Proton transfer continues by transferring H(6) to O(7) of
a third water molecule, finally producing an H
3
O
+
ion.
83
2. When two CO
2
molecules come close to each other as in Fig. 5.5(b), the
tendency toward C-C bond formation assists the ionization of a neighboring water
molecule. Here, by the transfer of H(2) to O(1), an OH
-
ion forms.
3. In the third pathway, shown in Fig. 5.5(c), transfer of H(2) to C(1) of a CO
2
molecule sandwiched between two water molecules generates an OH
-
ion.
4. Finally, in the forth pathway, shown in Fig. 5.5(d), transfer of H(2) to O(1) of a
CO
2
molecule leaves an OH
-
ion.
Fig. 5.5: Four dominant pathways to splitting of H
2
O molecules observed in MD
simulations. Each row displays one pathway, where the reactants and intermediate states
shown respectively in the left and right panels. Carbon, hydrogen and oxygen atoms are
shown respectively in green, white and red.
84
We should note that the pathways depicted in Fig. 5.5 are often part of long
chains of reactions. However, only key segments of such chains are shown for clarity of
presentation. For example, the third pathway shown in Fig. 5.5(c) is part of a chain, in
which a formic acid molecule (HCOOH) forms; see Fig. 5.6. Successive proton transfers
through these long chains are expected to prevent the recombination of H
+
and OH
-
ions,
which is consistent with Fig. 5.2(b). Namely, the number of ions is not decreasing at least
for last 1 ps of the simulation.
Fig. 5.6: A pathway for formation of a formic acid molecule. Each row displays a
reaction, left and right panels being reactants and products, respectively. These two
reactions take place after the reaction shown in Fig. 5.5(c). The formic acid molecule is
seen in panel (d). Transfers in panel (a) are: O (2) and H (3) to C (1), H (4) to O (5), H (6)
to O (7), H (8) to O (9). Moreover, unstable bond between H (10) and H (11) breaks. In
panel (c), H (2) transfers to O (1) and bond between O (1) and C (3) breaks.
85
We should also note that these pathways are not equally probable. For example,
the third pathway is by an order of magnitude less frequent. Moreover, we find the first
and the fourth pathways to have roughly the same likelihood. Although the first pathway
generates positive ions, the number of positive ions in the primary system is
approximately the same as that in the reference system (see Fig. 5.2(b)) indicating that
positive ions produced from the first pathway enter into new reactions. Whenever a π-
bond between the carbon atom and an oxygen atom of a CO
2
molecule breaks and the
carbon/oxygen atom forms a bond with OH/H from a water molecule (the first/fourth
pathway), another bond between the oxygen/carbon atom and a positive/negative ion
forms later. Therefore, the first and the forth pathways cannot be the reason for the excess
of negative ions in the primary system. We believe the excess of negative ions is a
manifestation of formation of C-C bonds as can be seen from the third pathway. The
average charge per ion for negative and positive ions is almost the same. Since there are
more negative ions than positive ions (in the primary system), the total charge of negative
ions (magnitude) is more than the total charge of positive ions, and the difference is
compensated by positive charges on molecules containing carbon atoms.
Fig. 5.7(a) shows the number of CO
2
molecules as a function of time in the
primary system. When the shock front is in the center of the nanobubble at t = 1.4 ps,
H
2
O-CO
2
reactions start, and consequently the number of CO
2
molecules start to decrease,
while the numbers of positive and negative ions start to increase (see Fig. 5.2(b)). Shortly
after the complete collapse of the nanobubble at t =2.3 ps (as the water hammer shock
86
spreads and local high pressure starts to decrease), some unstable products decay and lead
to a rise in the number of CO
2
molecules. However, most reactions are not reversible, at
least under the shock condition considered here. Among the products of H
2
O-CO
2
reactions are a small number of oxygen molecules, 8 of which are stable within the
duration of our simulation. Fig. 5.7(b) shows these oxygen molecules along with carbon
atoms.
Fig. 5.7: (a) Number of CO2 molecules as a function of time in the primary system. At t
= 1.4 ps, H2O-CO2 reactions start, and some unstable products start decaying at t = 2.3
ps. (b) Front view (the direction of shock propagation is perpendicular to the figure and is
outward) of carbon atoms and oxygen molecules at the end of simulation at t = 3.5 ps in
the primary system. Green spheres depict carbon atoms, and oxygen molecules are shown
in red.
We find that both H
2
O and CO
2
molecules act as sources of the oxygen atoms
constituting these oxygen molecules. However, the absence of any oxygen molecule in
the reference system implies that the presence of CO
2
molecules is necessary for the
formation of oxygen molecules from water molecules. The snapshots in Fig. 5.8 show the
formation of an oxygen molecule from a hydrogen peroxide (H
2
O
2
) molecule observed in
87
our simulation. Fig. 5.8(a) shows an H
2
O
2
molecule in the vicinity of a 2-carbon
molecule (the H
2
O
2
molecule consists of atoms originating from water molecules). Here,
an O atom (1) and C(2) are very close to each other. First, O(3) from the 2-carbon
molecule and H(4) from the H
2
O
2
molecule separate from their molecules and form a
hydroxyl ion; see Fig. 5.8(b). Subsequently, the second hydrogen atom H(5) of the H
2
O
2
molecule hops to a neighboring hydroxyl ion, forming a water molecule and leaving an
oxygen molecule, as shown in Fig. 5.8(c). We note that the presence of the 2-cabon
molecule surrounded by closely separated water molecules in the vicinity of the H
2
O
2
molecule is the key for triggering the reactions. Since this environment is created by the
high pressure associated with shock-induced collapse of the nanobubble, one may refer to
the reactions reported here by mechano-chemical reactions. Fig. 5.9 shows the formation
of another oxygen molecule, where one of the oxygen atoms initially belonged to a CO
2
molecule.
Fig. 5.8: Formation of an oxygen molecule from an H
2
O
2
molecule. Carbon, hydrogen
and oxygen atoms are shown in green, white, and red, respectively. Narrow sticks
between oxygen and hydrogen atoms indicate hydrogen bonds.
88
Fig. 5.9: Molecules involved in formation of an oxygen molecule. Carbon, hydrogen, and
oxygen atoms are shown in green, white, and red, respectively. Narrow bars between
oxygen and hydrogen atoms indicate hydrogen bonds. Here, a 2-carbon molecule is
surrounded by 2 water molecules and a water dimer. Oxygen atoms marked by 1 and 2
separate from their molecules and form an oxygen molecule.
89
6 Conclusions
The structure of water under shock, the dynamics of shock-induced collapse of
nanobubbles in water, and two promising applications of shock-induced nanobubble
collapse have been investigated. Firstly, the structure of water under shock has been
studied by a series of reactive force field MD simulations, which account for the
formation and dissociation of chemical bonds. The observed Hugoniot curve is in good
agreement with experimental results. We find molecular clusters with ice-VII-like
structure at particle velocity of 1 km/s. This is well supported by experimental
observations of ice VII for particle velocities between 0.75 and 2 km/s. Secondly, the
dynamics of shock-induced collapse of nanobubbles in water has been studied in detail.
When a shock wave strikes a nanobubble, the nanobubble starts shrinking. During
shrinkage, we observe a focused nanojet whose length scales linearly with the
nanobubble radius. This scaling relation has also been observed experimentally in shock-
induced collapse of millimeter-to-micron sized bubbles. We have examined the effects of
shock amplitude and the initial nanobubble size on the dynamics of nanobubble shrinkage
and collapse. The size of the nanojet increases with an increase in the nanobubble
diameter and shock amplitude. However the effect of nanobubble size is more
pronounced. Upon impact of the nanojet on the distal side of the nanobubble, a secondary
water hammer shock generates, which propagates spherically backward with a speed
higher than that of the primary shock. At the onset of its formation, the water hammer
shock has a pressure larger than the primary shock, but its pressure drops rapidly as it
90
spreads. Moreover, we observe a substantial number of positive and negative ions when
the nanojet hits the distal side of the nanobubble and water hammer shock generates.
As the first potential application of shock-induced collapse of nanobubbles,
poration in a lipid bilayer located in front of a collapsing nanobubble has been explored.
Our multimillion-atom classical MD simulations reveal that the nanojet impact
significantly deforms and thins the lipid bilayer and water molecules in the nanojet form
a spreading flow pattern after impact. Deformation and thinning of the lipid bilayer
combined with large pressure gradients across and spreading flow around the bilayer
create transient nanochannels through which water molecules translocate across the
bilayer.
Finally, we have studied mechano-chemical pathways to splitting of H
2
O and CO
2
molecules triggered by shock-induced nanobubble collapse. Without CO
2
molecules
inside the nanobubble, equal numbers of positive and negative ions form as a result of the
high pressure caused by the water hammer shock wave. In the presence of CO
2
molecules
inside the nanobubble, we observe an excess of negative ions. The difference arises from
H
2
O-CO
2
interactions, which start after the formation of the nanojet and before the
generation of the water hammer shock wave. A few oxygen and formic acid molecules
are among the products of the reactions. Proximity of molecules and enhanced hydrogen-
bond networks are found to play a crucial role in the splitting of H
2
O and CO
2
molecules
by shock-induced collapse of nanobubble. Dominant pathways, which lead to splitting of
water molecules, have been identified.
91
References
Abanades, S., A. Legal, et al. (2010). "Investigation of reactive cerium-based oxides for
H-2 production by thermochemical two-step water-splitting." Journal of Materials
Science 45(15): 4163-4173.
Ajaev, V. S. and G. M. Homsy (2006). "Modeling shapes and dynamics of confined
bubbles." Annual Review of Fluid Mechanics 38: 277-307.
Albon, N. and J. M. Sturtevant (1978). "Nature of Gel to Liquid-Crystal Transition of
Synthetic Phosphatidylcholines." Proceedings of the National Academy of
Sciences 75(5): 2258-2260.
Alder, B. J. and T. E. Wainwright (1957). "Phase transition for a hard sphere system."
Journal of Chemical Physics 27: 1208-1209.
Allen, M. P. and D. J. Tildesley (1987). Computer Simulation of Liquids. New York,
Oxford University Press.
Berendsen, H. J. C., J. P. M. Postma, et al. (1981). Intermolecular Forces.
Berendsen, H. J. C., J. P. M. Postma, et al. (1984). "Molecular-Dynamics with Coupling
to an External Bath." Journal of Chemical Physics 81(8): 3684-3690.
Berger, O., O. Edholm, et al. (1997). "Molecular dynamics simulations of a fluid bilayer
of dipalmitoylphosphatidylcholine at full hydration, constant pressure, and
constant temperature." Biophysical Journal 72(5): 2002-2013.
Berweger, C. D., W. F. Vangunsteren, et al. (1995). "Force-Field Parametrization by
Weak-Coupling - Reengineering Spc Water." Chemical Physics Letters 232(5-6):
429-436.
Borkent, B. M., S. M. Dammer, et al. (2007). "Superstability of surface nanobubbles."
Physical Review Letters 98(20): 204502.
Brennen, C. E. (1995). Cavitation and Bubble Dynamics. New York, Oxford University
Press.
Brujan, E. A., G. S. Keen, et al. (2002). "The final stage of the collapse of a cavitation
bubble close to a rigid boundary." Physics of Fluids 14(1): 85-92.
Cheung, K. S. and S. Yip (1991). "Atomic-Level Stress in an Inhomogeneous System."
Journal of Applied Physics 70(10): 5688-5690.
Chueh, W. C., C. Falter, et al. (2010). "High-Flux Solar-Driven Thermochemical
Dissociation of CO2 and H2O Using Nonstoichiometric Ceria." Science
330(6012): 1797-1801.
Chueh, W. C. and S. M. Haile (2009). "Ceria as a Thermochemical Reaction Medium for
Selectively Generating Syngas or Methane from H2O and CO2." Chemsuschem
2(8): 735-739.
92
Chueh, W. C. and S. M. Haile (2010). "A thermochemical study of ceria: exploiting an
old material for new modes of energy conversion and CO2 mitigation."
Philosophical Transactions of the Royal Society a-Mathematical Physical and
Engineering Sciences 368(1923): 3269-3294.
Didenko, Y. T. and K. S. Suslick (2002). "The energy efficiency of formation of photons,
radicals and ions during single-bubble cavitation." Nature 418(6896): 394-397.
Fisher, M. E. (1964). "The Free Energy of a Macroscopic System." Archive for Rational
Mechanics and Analysis 17(5): 377-410.
Flannigan, D. J. and K. S. Suslick (2005). "Plasma formation and temperature
measurement during single-bubble cavitation." Nature 434(7029): 52-55.
Frenkel, D. and B. Smit (2002). Understanding Molecular Simulation: From Algorithms
to Applications. San Diego, Academic Press.
Gaitan, D. F. and R. G. Holt (1999). "Experimental observations of bubble response and
light intensity near the threshold for single bubble sonoluminescence in an air-
water system." Physical Review E 59(5): 5495-5502.
Geissler, P. L., C. Dellago, et al. (2001). "Autoionization in liquid water." Science
291(5511): 2121-2124.
Golberg, A. and B. Rubinsky (2010). "A statistical model for multidimensional
irreversible electroporation cell death in tissue." Biomedical Engineering Online
9: 13.
Greengard, L. and V. Rokhlin (1997). "A fast algorithm for particle simulations
(Reprinted from the Journal of Computational Physics, vol 73, pg 325-348,
1987)." Journal of Computational Physics 135(2): 280-292.
Gstoehl, D., A. Brambilla, et al. (2008). "A quenching apparatus for the gaseous products
of the solar thermal dissociation of ZnO." Journal of Materials Science 43(14):
4729-4736.
Hansen, J. P. and I. R. McDonald (1986). Theory of Simple Liquids. New York,
Academic Press.
Hao, Y. and A. Prosperetti (1999). "The effect of viscosity on the spherical stability of
oscillating gas bubbles." Physics of Fluids 11(6): 1309-1317.
Hess, B., H. Bekker, et al. (1997). "LINCS: A linear constraint solver for molecular
simulations." Journal of Computational Chemistry 18(12): 1463-1472.
Hess, B., C. Kutzner, et al. (2008). "GROMACS 4: Algorithms for highly efficient, load-
balanced, and scalable molecular simulation." Journal of Chemical Theory and
Computation 4(3): 435-447.
Hilgenfeldt, S., D. Lohse, et al. (1996). "Phase diagrams for sonoluminescing bubbles."
Physics of Fluids 8(11): 2808-2826.
Hockney, R. W. and J. W. Eastwood (1981). Computer Simulation Using Particles. New
York, McGraw-Hill.
Holian, B. L. (1995). "Atomistic Computer-Simulations of Shock-Waves." Shock Waves
5(3): 149-157.
Holt, R. G. and D. F. Gaitan (1996). "Observation of stability boundaries in the parameter
space of single bubble sonoluminescence." Physical Review Letters 77(18): 3791-
3794.
93
Hong, K. S., H. F. Xu, et al. (2010). "Direct water splitting through vibrating
piezoelectric microfibers in water." Journal of Physical Chemistry Letters 1(6):
997-1002.
Hoover, W. G. (1985). "Canonical Dynamics - Equilibrium Phase-Space Distributions."
Physical Review A 31(3): 1695-1697.
Hoover, W. G. (1986). "Constant-Pressure Equations of Motion." Physical Review A
34(3): 2499-2500.
Janssens, G. O. A., B. G. Baekelandt, et al. (1995). "Comparison of Cluster and Infinite
Crystal Calculations on Zeolites with the Electronegativity Equalization Method
(Eem)." Journal of Physical Chemistry 99(10): 3251-3258.
Johnsen, E. and T. Colonius (2008). "Shock-induced collapse of a gas bubble in
shockwave lithotripsy." Journal of the Acoustical Society of America 124(4):
2011-2020.
Johnston, H. S. and C. Parr (1963). "Activation Energies from Bond Energies .1.
Hydrogen Transfer Reactions." Journal of the American Chemical Society 85(17):
2544-&.
Jorgensen, J. D. and T. G. Worlton (1985). "Disordered Structure of D2o Ice Vii from
Insitu Neutron Powder Diffraction." Journal of Chemical Physics 83(1): 329-333.
Kaneko, H., T. Miura, et al. (2007). "Reactive ceramics of CeO2-MOx (M = Mn, Fe, Ni,
Cu) for H-2 generation by two-step water splitting using concentrated solar
thermal energy." Energy 32(5): 656-663.
Ketterling, J. A. and R. E. Apfel (2000). "Extensive experimental mapping of
sonoluminescence parameter space." Physical Review E 61(4): 3832-3837.
Ketterling, J. A. and R. E. Apfel (2000). "Shape and extinction thresholds in
sonoluminescence parameter space." Journal of the Acoustical Society of
America 107(3): L13-L18.
Klotz, S., L. E. Bove, et al. (2009). "The preparation and structure of salty ice VII under
pressure." Nature Materials 8(5): 405-409.
Knuth, D. (1973). The Art of Computer Programming. Reading, MA, Addison-Wesley.
Kodama, T. and N. Gokon (2007). "Thermochernical cycles for high-temperature solar
hydrogen production." Chemical Reviews 107(10): 4048-4077.
Kodama, T. and K. Takayama (1998). "Dynamic behavior of bubbles during
extracorporeal shock-wave lithotripsy." Ultrasound in Medicine and Biology
24(5): 723-38.
Kodama, T. and Y. Tomita (2000). "Cavitation bubble behavior and bubble-shock wave
interaction near a gelatin surface as a study of in vivo bubble dynamics." Applied
Physics B-Lasers and Optics 70(1): 139-149.
Kudo, N., K. Okada, et al. (2009). "Sonoporation by Single-Shot Pulsed Ultrasound with
Microbubbles Adjacent to Cells." Biophysical Journal 96(12): 4866-4876.
Kuijpers, M. W. A., D. van Eck, et al. (2002). "Cavitation-induced reactions in high-
pressure carbon dioxide." Science 298(5600): 1969-1971.
Lebowitz, J. L., J. K. Percus, et al. (1967). "Ensemble Dependence of Fluctuations with
Application to Machine Computations." Physical Review 153(1): 250-&.
94
Lewis, N. S. and D. G. Nocera (2006). "Powering the planet: Chemical challenges in
solar energy utilization." Proceedings of the National Academy of Sciences of the
United States of America 103(43): 15729-15735.
Lin, H., B. D. Storey, et al. (2002). "Rayleigh-Taylor instability of violently collapsing
bubbles." Physics of Fluids 14(8): 2925-2928.
Martyna, G. J., M. L. Klein, et al. (1992). "Nose-Hoover Chains - the Canonical
Ensemble Via Continuous Dynamics." Journal of Chemical Physics 97(4): 2635-
2643.
McQuarrie, D. A. (1976). Statistical Mechanics. New York, Harper and Row.
Metropolis, N., A. W. Rosenbluth, et al. (1953). "Equation of State Calculations by Fast
Computing Machines." Journal of Chemical Physics 21(6): 1087-1092.
Miller, J. E., M. D. Allendorf, et al. (2008). "Metal oxide composites and structures for
ultra-high temperature solar thermochemical cycles." Journal of Materials Science
43(14): 4714-4728.
Mitragotri, S. (2005). "Innovation - Healing sound: the use of ultrasound in drug delivery
and other therapeutic applications." Nature Reviews Drug Discovery 4(3): 255-
260.
Mortier, W. J., S. K. Ghosh, et al. (1986). "Electronegativity Equalization Method for the
Calculation of Atomic Charges in Molecules." Journal of the American Chemical
Society 108(15): 4315-4320.
Mortier, W. J., K. Vangenechten, et al. (1985). "Electronegativity Equalization -
Application and Parametrization." Journal of the American Chemical Society
107(4): 829-835.
Moseler, M. and U. Landman (2000). "Formation, stability, and breakup of nanojets."
Science 289(5482): 1165-1169.
Nakano, A., R. K. Kalia, et al. (2007). "A divide-and-conquer/cellular-decomposition
framework for million-to-billion atom simulations of chemical reactions."
Computational Materials Science 38(4): 642-652.
Newman, C. M. H. and T. Bettinger (2007). "Gene therapy progress and prospects:
Ultrasound for gene transfer." Gene Therapy 14(6): 465-475.
Nicolas, J. J., K. E. Gubbins, et al. (1979). "Equation of State for the Lennard-Jones
Fluid." Molecular Physics 37(5): 1429-1454.
Nomura, K. I., R. K. Kalia, et al. (2008). "A scalable parallel algorithm for large-scale
reactive force-field molecular dynamics simulations." Computer Physics
Communications 178(2): 73-87.
Nomura, K. I., R. K. Kalia, et al. (2007). "Dynamic transition in the structure of an
energetic crystal during chemical reactions at shock front prior to detonation."
Physical Review Letters 99(14): 046304-046308.
Nose, S. (1984). "A Molecular-Dynamics Method for Simulations in the Canonical
Ensemble." Molecular Physics 52(2): 255-268.
Nose, S. (1984). "A Unified Formulation of the Constant Temperature Molecular-
Dynamics Methods." Journal of Chemical Physics 81(1): 511-519.
Ohl, C. D., M. Arora, et al. (2006). "Surface cleaning from laser-induced cavitation
bubbles." Applied Physics Letters 89(7): 074102.
95
Ohl, C. D., M. Arora, et al. (2006). "Sonoporation from jetting cavitation bubbles."
Biophysical Journal 91(11): 4285-4295.
Ohl, C. D. and R. Ikink (2003). "Shock-wave-induced jetting of micron-size bubbles."
Physical Review Letters 90(21): 214502.
Ohl, S. W., E. Klaseboer, et al. (2009). "The dynamics of a non-equilibrium bubble near
bio-materials." Physics in Medicine and Biology 54(20): 6313-6336.
Parr, R. G. and R. G. Pearson (1983). "Absolute hardness: companion parameter to
absolute electronegativity." Journal of the American Chemical Society 105(26):
7512-7516.
Perkins, C. and A. W. Weimer (2009). "Solar-Thermal Production of Renewable
Hydrogen." Aiche Journal 55(2): 286-293.
Philipp, A. and W. Lauterborn (1998). "Cavitation erosion by single laser-produced
bubbles." Journal of Fluid Mechanics 361: 75-116.
Plesset, M. S. and R. B. Chapman (1971). "Collapse of an Initially Spherical Vapour
Cavity in Neighbourhood of a Solid Boundary." Journal of Fluid Mechanics
47(May31): 283-&.
Powles, J. G., W. A. B. Evans, et al. (1982). "Non-Destructive Molecular-Dynamics
Simulation of the Chemical-Potential of a Fluid." Molecular Physics 46(6): 1347-
1370.
Prosperetti, A. (1977). "Viscous Effects on Perturbed Spherical Flows." Quarterly of
Applied Mathematics 34(4): 339-352.
Quentrec, B. and C. Brot (1973). "New Method for Searching for Neighbors in Molecular
Dynamics Computations." Journal of Computational Physics 13(3): 430-432.
Rappe, A. K. and W. A. Goddard (1991). "Charge Equilibration for Molecular-Dynamics
Simulations." Journal of Physical Chemistry 95(8): 3358-3363.
Roeb, M. and H. Muller-Steinhagen (2010). "Concentrating on solar electricity and
fuels." Science 329(5993): 773-774.
Root, D. M., C. R. Landis, et al. (1993). "Valence Bond Concepts Applied to the
Molecular Mechanics Description of Molecular Shapes .1. Application to
Nonhypervalent Molecules of the P-Block." Journal of the American Chemical
Society 115(10): 4201-4209.
Rosenthal, I., J. Z. Sostaric, et al. (2004). "Sonodynamic therapy - a review of the
synergistic effects of drugs and ultrasound." Ultrasonics Sonochemistry 11(6):
349-363.
Rybakov, A. P. and I. A. Rybakov (1995). "Polymorphism of Shocked Water." European
Journal of Mechanics B-Fluids 14(3): 323-332.
Sanderson, R. T. (1976). Chemical Bond and Bond Energies. New York, Academic Press.
Sanderson, R. T. (1983). Polar Covalence. New York, Academic Press.
Sankin, G. N. and P. Zhong (2006). "Interaction between shock wave and single inertial
bubbles near an elastic boundary." Physical Review E 74(4): 046304-046308.
Schunk, L. O., W. Lipinski, et al. (2009). "Ablative Heat Transfer in a Shrinking Packed-
Bed of ZnO Undergoing Solar Thermal Dissociation." Aiche Journal 55(7): 1659-
1666.
96
Schunk, L. O. and A. Steinfeld (2009). "Kinetics of the Thermal Dissociation of ZnO
Exposed to Concentrated Solar Irradiation Using a Solar-Driven
Thermogravimeter in the 1800-2100 K Range." Aiche Journal 55(6): 1497-1504.
Steinfeld, A. (2005). "Solar thermochemical production of hydrogen - a review." Solar
Energy 78(5): 603-615.
Stoddard, S. D. and J. Ford (1973). "Numerical Experiments on Stochastic Behavior of a
Lennard-Jones Gas System." Physical Review A 8(3): 1504-1512.
Strachan, A., A. C. T. van Duin, et al. (2003). "Shock waves in high-energy materials:
The initial chemical events in nitramine RDX." Physical Review Letters 91(9):
098301.
Sun, Z. H., X. X. Wang, et al. (2006). "On stress calculations in atomistic simulations."
Modelling and Simulation in Materials Science and Engineering 14(3): 423-431.
Tieleman, D. P. and H. J. C. Berendsen (1996). "Molecular dynamics simulations of a
fully hydrated dipalmitoyl phosphatidylcholine bilayer with different macroscopic
boundary conditions and parameters." Journal of Chemical Physics 105(11):
4871-4880.
Tomita, Y. and T. Kodama (2003). "Interaction of laser-induced cavitation bubbles with
composite surfaces." Journal of Applied Physics 94(5): 2809-2816.
Tomita, Y. and A. Shima (1986). "Mechanisms of Impulsive Pressure Generation and
Damage Pit Formation by Bubble Collapse." Journal of Fluid Mechanics 169:
535-564.
Tsai, D. H. (1979). "Virial Theorem and Stress Calculation in Molecular-Dynamics."
Journal of Chemical Physics 70(3): 1375-1382.
Tuckerman, M., B. J. Berne, et al. (1992). "Reversible Multiple Time Scale Molecular-
Dynamics." Journal of Chemical Physics 97(3): 1990-2001.
Tufaile, A. and J. C. Sartorelli (2002). "Bubble and spherical air shell formation
dynamics." Physical Review E 66(5): 056204.
Turner, J. A. (2004). "Sustainable hydrogen production." Science 305(5686): 972-974.
van Duin, A. C. T., S. Dasgupta, et al. (2001). "ReaxFF: A reactive force field for
hydrocarbons." Journal of Physical Chemistry A 105(41): 9396-9409.
van Duin, A. C. T., A. Strachan, et al. (2003). "ReaxFF(SiO) reactive force field for
silicon and silicon oxide systems." Journal of Physical Chemistry A 107(19):
3803-3811.
Vangenechten, K. A., W. J. Mortier, et al. (1987). "Intrinsic Framework Electronegativity
- a Novel Concept in Solid-State Chemistry." Journal of Chemical Physics 86(9):
5063-5071.
Vernier, P. T., Y. H. Sun, et al. (2006). "Nanoelectropulse-driven membrane perturbation
and small molecule permeabilization." Bmc Cell Biology 7: 37.
Warsi, Z. U. A. (1998). Fluid Dynamics: Theoretical and Computational Approaches,
CRC Press.
Wood, W. W. and J. D. Jacobson (1957). "Preliminary results from a recalculation of the
Monte Carlo equation of state of hard-spheres." Journal of Chemical Physics 27:
1207-1208.
97
Yamanoi, I. and M. Tamagawa (2006). "Deformation analysis of bubble near curved
elastic wall for developing shock wave DDS." Jsme International Journal Series
B-Fluids and Thermal Engineering 49(3): 755-760.
Zhou, Y., J. M. Cui, et al. (2008). "Dynamics of sonoporation correlated with acoustic
cavitation activities." Biophysical Journal 94(7): L51-L53.
Zwaan, E., S. Le Gac, et al. (2007). "Controlled cavitation in microfluidic systems."
Physical Review Letters 98(25): 254501.
Abstract (if available)
Abstract
The shock-induced collapse of nanobubbles in water is investigated using molecular dynamics simulations based on a reactive force field. Monitoring the collapse of a cavitation nanobubble, we observe a focused nanojet at the onset of bubble shrinkage and a water hammer shock wave upon bubble collapse. The nanojet length scales linearly with the nanobubble radius, as observed in experiments on micron-to-millimeter size bubbles. The shock induces dramatic structural changes, including an ice-VII-like structural motif at a particle velocity of approximately 1 km/s. The incipient ice VII formation and the calculated Hugoniot curve are in good agreement with experimental results. Moreover, a substantial number of positive and negative ions appear when the nanojet hits the distal side of the nanobubble and the water hammer shock forms. Furthermore, two promising applications of shock-induced nanobubble collapse have been explored. Our simulations of poration in lipid bilayers due to shock-induced collapse of nanobubbles reveal penetration of nanojets into lipid bilayers. The nanojet impact generates shear flow of water on bilayer leaflets and pressure gradients across them, which transiently enhance the bilayer permeability by creating nanopores through which water molecules translocate across the bilayer. The effects of nanobubble size and temperature on the porosity of lipid bilayers are examined. Finally, the shock-induced collapse of CO2-filled nanobubbles in water is investigated. The energetic nanojet and high-pressure water hammer shock formed during and after collapse of the nanobubble trigger mechano-chemical H2O-CO2 reactions, some of which lead to splitting of water molecules. The dominant pathways through which splitting of water molecules occur are identified.
Linked assets
University of Southern California Dissertations and Theses
Conceptually similar
PDF
Shock-induced poration, cholesterol flip-flop and small interfering RNA transfection in a phospholipid membrane: multimillion atom, microsecond molecular dynamics simulations
PDF
Theoretical studies of lipid bilayer electroporation using molecular dynamics simulations
PDF
Molecular simulations of water and monovalent ion dynamics in the electroporation of phospholipid bilayers
PDF
A measurement of the dielectric properties of water nanoclusters
PDF
Modeling nanodevices: from semiconductor heterostructures to Josephson junction arrays
PDF
Modeling graphene: magnetic, transport and optical properties
PDF
Towards optimized dynamical error control and algorithms for quantum information processing
PDF
Converging shocks in water and material effects
PDF
An experimental investigation by optical methods of the physics and chemistry of transient plasma ignition
PDF
Microbially synthesized chalcogenide nanostructures: characterization and potential applications
PDF
Nanomaterials under extreme environments: a study of structural and dynamic properties using reactive molecular dynamics simulations
PDF
One-dimensional nanomaterials: synthesis and applications
PDF
Shock sensitivity of energetic material and nanometric damage mechanisms in silica glass
PDF
Electronic correlation effects in multi-band systems
PDF
Entanglement parity effects in quantum spin chains
PDF
Modeling the hot-dense plasma of the solar interior in and out of thermal equilibrium
PDF
Bacterial nanowires of Shewanella oneidensis MR-1: electron transport mechanism, composition, and role of multiheme cytochromes
PDF
The physics of membrane protein polyhedra
PDF
Quantum computation by transport: development and potential implementations
PDF
Quantum phase transitions in disordered antiferromagnets
Asset Metadata
Creator
Vedadi, Mohammad Hossein
(author)
Core Title
Shock-induced nanobubble collapse and its applications
School
College of Letters, Arts and Sciences
Degree
Doctor of Philosophy
Degree Program
Physics
Publication Date
03/04/2013
Defense Date
02/07/2013
Publisher
University of Southern California
(original),
University of Southern California. Libraries
(digital)
Tag
bubble collapse,molecular dynamics,OAI-PMH Harvest
Language
English
Contributor
Electronically uploaded by the author
(provenance)
Advisor
Haas, Stephan (
committee chair
), Dappen, Werner (
committee member
), El-Naggar, Mohammed (
committee member
), Kresin, Vitaly (
committee member
), Ma, Jin (
committee member
)
Creator Email
mohammad.vedadi@gmail.com,vedadi@usc.edu
Permanent Link (DOI)
https://doi.org/10.25549/usctheses-c3-222938
Unique identifier
UC11292763
Identifier
usctheses-c3-222938 (legacy record id)
Legacy Identifier
etd-VedadiMoha-1456.pdf
Dmrecord
222938
Document Type
Dissertation
Rights
Vedadi, Mohammad Hossein
Type
texts
Source
University of Southern California
(contributing entity),
University of Southern California Dissertations and Theses
(collection)
Access Conditions
The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the a...
Repository Name
University of Southern California Digital Library
Repository Location
USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA
Tags
bubble collapse
molecular dynamics