ATOMIC-SCALE DYNAMIC PROCESSES

IN THE BRITTLE FRACTURE OF SILICA

By

THOMAS P. SWILER

A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL

OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT

OF THE REQUIREMENTS FOR THE DEGREE OF

DOCTOR OF PHILOSOPHY

UNIVERSITY OF FLORIDA

Copyright 1994

by

Thomas P. Swiler

This dissertation is dedicated to my parents, Thomas Russell Swiler and Elaine

Palda Swiler, who were always forgiving of my early "fracture experiments." Until we

can fulfill our destructive tendencies through simulation, we must rely on forgiveness.

ACKNOWLEDGEMENTS

I would like to acknowledge the support of my advisor, Joseph H. Simmons,

throughout this work. He has continually supported a molecular dynamics effort even

though it has not been supported by external sources. The opportunity to explore the

implications of this simulation work without having to provide specific answers to

questions posed by sponsors has made work in this area interesting and instructive.

I would like to thank John J. Mecholsky, Jr. for enlightening me to current

research in brittle fracture. His observations and interpretations of fractal fracture

surfaces have provided me with the background to make the link between computer

simulated fracture and experimentally observed fracture.

I would like to thank Adrian C. Wright for performing comparisons of our pair

correlation function data derived from molecular dynamic simulations to his radial

distribution data derived from neutron diffraction experiments.

Finally, I would like to thank my family and friends who have provided support

for me during the up and downs of writing this dissertation.

TABLE OF CONTENTS

ACKNOWLEDGEMENTS .............. ................... .. iv

LIST OF TABLES ..................................... ix

LIST OF FIGURES .................. .................... x

KEY TO SYMBOLS ............ .............. ........... xv

Abstract .............. .................. ............ xxiii

CHAPTERS

1 INTRODUCTION ............... .......... ........ 1

W hat is Fracture .......................... ......... 1

Fracture as a Transformation .............. ... ............ 2

Molecular Dynamics ................... ............ 4

Statement of Problem ............. .................... 7

2 BACKGROUND ............. ............ .......... 8

Brittle Fracture ............... ................... 9

Characteristics of Brittle Fracture ........... ..... ..... 9

The Ultimate Strength of Materials ........... ..... .... 9

The Effect of Flaws ........... ...... ......... . 11

Regimes of Brittle Fracture ............. .... ....... 12

Environmental Effects ................ . .......... 12

Crack Propagation Velocities ............. .... ....... 14

Reading the Fracture Surface .......... ...... ..... 14

Silica ............................... ........ 16

Bonding in Silica ...................... ........ 16

Structure of Silica ................... ...... .... 16

Phases of Silica .............................. 18

The Crystal Structure of Cristobalite ............... ... 19

The Ring Structure of Silica ............ ..... ....... 20

Elastic Properties of Vitreous Silica . . . . . . ..... 23

The Fracture Strength of Silica . . . . . . . . .. 24

Stress Corrosion of Silica . . . . . ..... ........ 25

Molecular Dynamics ................... ....... ... . . 28

Survey of Applications .................... .... .. . 28

Survey of Potential Functions . . . . . . . ..... . ... 31

Potentials for Silica . . . . . ........ .. .. . . 34

Applications of MD in Models of Silica . . . . . ..... 37

Comparisons of Simulated Structure to Experiment . . . ... 39

MD Simulations of Fracture . . . . . ... . 40

Chaos in Fracture: Determinism and Fractal Surfaces .......... . 43

Investigations of Chaotic Processes . . . . . . . ..... 44

Predictability of Brittle Fracture . . . . . . . . 46

Fractal Geometry ................... ........... 47

Applications of Fractal Geometry to Fracture . . . . .... 49

The Fractal-like Nature of Fracture Surfaces . . . . .... 51

Applicability of Fractal Analyses to MD Simulated Systems . . 52

3 PROCEDURE .................. ......... ......... 56

Using Molecular Dynamics to Model Fracture . . . . . . . . 56

Potential Functions ........... .... ... .......... 58

The Bom-Mayer-Huggins Potential . . . . . . . ..... 60

Truncation or Screening Functions . . . . . . . ..... 61

The Soules Potential Function . . . . . . . . .. 63

The Feuston-Garofalini Potential Function . . . . . ..... 65

Simulation Techniques ........ . . . . . . . . ...... 68

Integration Algorithms ... . . . . . . . ....... 68

Verlet algorithm ................... ...... 68

Fifth-order Gear predictor-corrector algorithm . . ... 69

The Time Step ........... ...... .. ......... 69

Limiting Particle Interactions . . . . . . . 70

Application of Periodic Boundary Conditions . . . . . ... 72

Handling Particle Interactions . . . . . . . . 74

Observing the System . . . . . ......... . . .. 76

System thermodynamics . . . . . . . . 76

System structure .................. ........ 77

Force Distributions . . . . . . . .. 79

Damped Motion Molecular Dynamics . . . . . . . ... 80

Assigning Initial Velocities . . . . . . . . .. 80

Analytical Techniques . . . . . . ....... . . . ..... 81

Graphic Representations . . . . . . . . .. 82

Free Volume Sphere Analysis . . . . . . . .. 82

Ring Statistics ...................... .......... 84

Structural Defects ................. ......... 86

Composition Profiles ........................... 86

Mechanical Properties .......................... 87

Program Details ................ . . . . . . . .... 87

Overview .................................. 87

The Structure of Computer Code . . . . . . ...... 88

The Choice of a Programming Language . . . . . ..... 92

The Simulation Program . . . . .......... ...... 95

The Graphic Analysis Program . . . . . . . .... . 96

Other Analysis Programs . . . . . ....... . . . 97

Experimental Techniques ................... . . . . 97

Initial Conditions ................... ........... 97

Application of Uniaxial Strain . . . . . . . ... . 99

Application of Shear Strain ........... . . . . 101

Experiments and Analyses . . . . . ........ . . . . ..... 103

Phases Simulated ................... . . . . 103

Strain Rates Simulated . . . . . . . . . ...... 103

Comparison of Potential Functions .. . . . . . . . 103

Thermal Effects ............................. 104

Studies of Determinism . . . . .......... ........ 104

Studies of Nonlinear Viscosity . . . . . . . 105

4 RESULTS AND DISCUSSION ................... ...... 106

The Models ............................. ........ 106

Static (Structural) Characteristics . . . . . . . 107

Densities and pressure tensors . . . . . . . ... 107

Low-order pair-correlation functions . . . . . .... 108

High-order pair correlation functions . . . . . ... 112

Discussion of higher order peaks . . . . . . ... 113

Bond angle distributions . . . . . . . . . . 116

Comparison of structural distributions to experiment . .. 117

Ring statistics ........ .................. 121

Volume distributions . . . . . . . . ... 124

Defects ......... ..................... 125

Dynamic Properties ............... ............. 126

Melting temperature . . . . . . ... ....... 126

Mechanical properties . . . . . . . . ... 129

Review of Models .................. ......... 131

Fracture by Application of Uniaxial Strain . . . . . . . . ... 132

Overview ................. ........ .......... 132

Low Cristobalite ................... ........... 133

Stress-strain curves . . . . ......... . . 134

Energetics of fracture . . . . . . . . . 142

Fracture surfaces . . . . . . . . 149

Vitreous Silica. ...... .. ........... 154

Stress-strain curves . . . . . . . 154

Energetics of fracture ....................... 160

Structural changes during strain-to-fracture . . . 164

Viscous Flow ........... .............. ........... 174

Determinism in Fracture ....................... ... .. .. 179

Reproducibility and Accuracy of Computer Simulations . . . 179

Variation in Anneal Times ..... . . . . . . . . 180

Round-off Errors ....... .. . . . . . . . . . 182

"Butterfly Effect" ........... ... ............ .. 183

Sensitivity on Simulation Artifacts . . . . . . ..... 185

5 ANALYSIS ......... ................... ........ .. 187

M odels of Silica .................. . .. ........... 187

Stress and Strain ............... ...... ............. 190

Structural Changes During Strain-to-Failure . . . ... . . . 190

Shear Flow .......... .. . .... . . .......... 193

Initiation of Fracture . . . . . .. ... . . . . ...... 194

Crack Propagation .......... . . . . . . . . . 198

Simulation of Stochastic Crack Growth . . . . . . . . .. 199

The Completion of Fracture . . . . . ........ . . . . 202

The Structure of the Fracture Surface . . . . ... . . . ... 203

When is a Bond Broken ..... ........ ......... .... 203

Surface Reconstruction . . . . ........ . . . . 204

Determinism in Fracture ...... . . . . . . . . . . 207

6 CONCLUSIONS AND FUTURE WORK .......... . . . .. . 209

REFERENCE LIST ................ ...... .......... .. 211

BIOGRAPHICAL SKETCH ........ . . . . . . . . . 219

LIST OF TABLES

Table

2-1 Characteristics of the Phases of Silica . . . . . . . . . 19

3-1 Parameters Used in the Soules Force and Potential Functions ...... . 63

3-2 Parameters used in the Feuston-Garofalini Potential Function . . ... 67

4-1 Coordination of species in simulations using various potentials and

methods of formation as found by the ring statistics program. ...... 125

4-2 Coordination of species in simulations of vitreous silica at 0.2 strain and

various strain rates as found by the ring statistics program. . . ... 173

4-3 Viscosities of vitreous silica as modeled by the Soules potential at various

temperatures and strain rates, in units of 1021 Pas. . . . . .... 175

LIST OF FIGURES

Figure

2-1. Schematic showing crack propagation velocities as a function of stress

intensity factor.. ................. . . . . . . . . ...... 13

2-2. The silica tetrahedron is the basic structural unit of silica structures, a)

One silica tetrahedron. b) Two tetrahedra linked by comer sharing. 0 is

the Si-O-Si bond angle.. ................... . . . ... 17

2-3. Chemical reaction proposed by Michalske and Freiman for

environmentally enhanced Si-O-Si bond rupture. The strained crack tip

bond (a) adsorbs a water molecule, (b) is cleaved by dissociative reaction,

and (c) is converted to surface silanol groups.. . . . . . . ... 26

3-1. Plot of the Soules potential function for various interactions. . . ... 64

3-2. Plot of the two-body portion of the Feuston-Garofalini potential function

for various interactions. .................. ... . . . . 66

3-3. The effect of periodic boundary conditions in a simulated system. Note

that the neighborhood of an atom may extend beyond the primary cell. 71

3-4. Schematic showing the application of periodic boundary conditions to a

system under shear strain. ........ . . . . . . . . 74

4-1. Comparison the first peaks of the Si-O pair correlation functions for

vitreous silica as modeled by the Soules and the F-G potential

functions . . . . .. . . .. .. .. . . . . . . . . . 109

4-2. Comparison the first peaks of the 0-0 pair correlation functions for

vitreous silica as modeled by the Soules and the F-G potential

functions ................... ............ ......... 110

4-3. Comparison the first peaks of the Si-Si pair correlation functions for

vitreous silica as modeled by the Soules and the F-G potential

functions ................... ......... ........... 111

4-4. Comparison the higher order peaks of the Si-O pair correlation functions

for vitreous silica as modeled by the Soules and the F-G potential

functions ................... ................... 113

4-5. Contribution of various pair interactions to the potential energy of vitreous

silica as modeled by the Soules potential function. . . . . . ... 114

4-6. Contribution of various pair interactions to the potential energy of vitreous

silica as modeled by the F-G potential function. . . . . . .... 115

4-7. Comparison of O-Si-O bond angle distributions for vitreous silica as

modeled by the Soules and the F-G potential functions. . . . . ... 117

4-8. Comparison of Si-O-Si bond angle distributions for vitreous silica as

modeled by the Soules and the F-G potential functions. . . . . ... 118

4-9. Comparison of radial distribution functions derived from a neutron

diffraction experiment (solid line) and a molecular dynamics simulation

using the Soules potential (dashed line). . . . . . . ...... ...... 119

4-10. Comparison of radial distribution functions derived from a neutron

diffraction experiment (solid line) and a molecular dynamics simulation

using the F-G potential (dashed line). . . . . . . . . . 120

4-11. Comparison of the ring size distributions for vitreous silica as modeled by

the Soules and the F-G potential functions and as formed by crystal melt

and random position techniques. . . . . . ..... ....... 122

4-12. Comparison of the free volume sphere size distributions for vitreous silica

as modeled by the Soules and the F-G potential functions and as formed

by various techniques. ........ . . . . . . . . . 124

4-13. Potential energy as a function of temperature for cristobalite as modeled

by the Soules potential. ......... . . . . . . . ....... 127

4-14. Potential energy as a function of temperature for cristobalite as modeled

by the F-G potential .................. .............. 128

4-15. Equilibrium elastic stress-strain curves for cristobalite and vitreous silica

as modeled by the Soules and the F-G potential functions. . . . ... 129

4-16. Comparison of stress-strain curves at various strain rates for low

cristobalite as modeled by the Soules potential. . . . . . . ... 134

4-17. Comparison of stress-strain curves at various strain rates for low

cristobalite as modeled by the F-G potential. . . . . . . .... 135

4-18. View down the C-axis of low cristobalite under an intermediate tensile

strain applied along the C-axis, as modeled by the Soules potential. In

this view and in all following views, oxygen atoms are represented by the

larger, lighter balls, and silicon atoms are represented by the smaller,

darker bals. ..................... ........ ........... 137

4-19. View down the C-axis of low cristobalite under an intermediate tensile

strain applied along the C-axis, as modeled by the F-G potential. . . 138

4-20. View down the C-axis of low cristobalite under a high tensile strain

applied along the C-axis, as modeled by the F-G potential. . . . ... 139

4-21. View down the C-axis of low cristobalite under a high tensile strain

applied along the C-axis, as modeled by the Soules potential. ...... .140

4-22. Potential energy as a function of strain for low cristobalite as modeled by

the Soules potential function under adiabatic strain at various strain

rates. ................... .............. ........ 144

4-23. Kinetic energy as a function of strain for low cristobalite as modeled by

the Soules potential function under adiabatic strain at various strain

rates. ................... ............. ......... 145

4-24. System temperature as a function of strain for low cristobalite as modeled

by the Soules potential function under adiabatic strain at various strain

rates. ................... ............. ......... 146

4-25. Schematic showing the flow of energy in a brittle material system during

strain-to-failure ........... . . . . . . . . . ..... 147

4-26. Potential energy as a function of strain for low cristobalite as modeled by

the Soules potential function under isothermal strain at various strain

rates. ................... ............ .......... 148

4-27. Normal view of a shallow slice of a surface created during the fracture of

low cristobalite as modeled by the Soules potential under isothermal

conditions at 0.05 ps' strain rate. The failure of the slice to fill the entire

area is a consequence of the roughness of the surface. . . . . ... 149

4-28. Side view of surfaces created during the fracture of low cristobalite as

modeled by the Soules potential under isothermal conditions at 0.05 ps'

strain rate .......... ....... ...... .............. 151

4-29. Side view of surfaces created during the fracture of low cristobalite as

modeled by the Soules potential under isothermal conditions at 0.5 ps'-

strain rate .................. ...... .............. 151

4-30. Relative abundance of oxygen and silicon atoms in slices along the length

of the cell perpendicular to the direction of strain. Values plotted are the

predominance of oxygen over silicon. . . . . . . . . .. 152

4-31. Close-up side view of surface created during the fracture of low

cristobalite as modeled by the Soules potential under isothermal conditions

at 0.5 ps' strain rate ................... .. ........... 153

4-32. Comparison of stress-strain curves at various strain rates for vitreous silica

as modeled by the Soules potential. . . . . . . . . 155

4-33. Comparison of stress-strain curves at various strain rates for vitreous silica

as modeled by the F-G potential. ..... . . . . . . 156

4-34. Distribution of forces on the atoms comprising a vitreous silica sample

under strain at 0.05 ps' strain rate. Various curves represent different

engineering strains. ................. . . . . . 158

4-35. Distribution of forces on the atoms comprising a vitreous silica sample

under strain at 5.0 ps' strain rate. Various curves represent different

engineering strains. ......... . . . . . . . . . . 159

4-36. Potential energy as a function of strain for vitreous silica as modeled by

the Soules potential function under adiabatic strain at various strain

rates. ................................... ..... 161

4-37. System temperature as a function of strain for vitreous silica as modeled

by the Soules potential function under adiabatic strain at various strain

rates. ................... ..................... .. 162

4-38. Potential energy as a function of strain for vitreous silica as modeled by

the Soules potential function under isothermal strain at various strain

rates. ................... .................... .. 163

4-39. Si-O pair correlations in vitreous silica. Comparisons between unstrained

and strained to 0.2 at various strain rates are made. . . . . .... 165

4-40. Si-Si pair correlations in vitreous silica. Comparisons between unstrained

and strained to 0.2 at various strain rates are made. . . . . .... 166

4-41. O-Si-O bond angle distributions in vitreous silica. Comparisons between

unstrained and strained to 0.2 at various strain rates are made. ...... 167

4-42. Si-O-Si bond angle distributions in vitreous silica. Comparisons between

unstrained and strained to 0.2 at various strain rates are made. . . . 168

4-43. Comparison of ring size distributions of vitreous silica strained to 0.2 at

different strain rates. ........... . . . . . . . . 169

4-44. Change in the abundances of rings of various sizes during strain to 0.2 in

vitreous silica at various strain rates. . . . . . . . 170

4-45. Free volume sphere histogram--the integral of the free volume sphere

distribution--for vitreous silica strained to 0.2 strain at various strain

rates. ................... ........... . ........ . 171

4-46. Shear stress as a function of time for a vitreous silica sample as modeled

by the Soules potential function at various temperatures. . . . ... 176

4-47. Shear stress as a function of time for a vitreous silica sample as modeled

by the Soules potential function at various strain rates. . . . . ... 177

4-48. Stress-strain curves for nearly identical systems after various anneal

times. ........... .... ......... ............ 180

4-49. Stress-strain curves for nearly identical systems differing by a loss of

precision resulting from truncation of system conditions at various

times. ................ . ........ ........ .. 183

4-50. Deviations from stress-strain curves resulting from a displacement of one

atom immediately before straining. . . . . . . . . 184

4-51. Deviations from stress-strain curves resulting from a displacement of one

atom before an anneal period prior to straining. . . . . . . ... 185

5-1. Pseudo fracture surface created by simple stochastic crack propagation

simulation ......... .... .................. 201

5-2. Schematic showing the fracture of a bond in vitreous silica. Series shows

unstrained bond, strained bond, and broken bond. . . . . . ... 205

KEY TO SYMBOLS

A 1. A parameter in the Born potential. 2. The cross sectional area covered

by a bond. 3. A parameter in one form of the Lennard-Jones potential.

4. A parameter in the two-body portion of the Stillinger-Weber potential.

A,,, A parameter in the two-body part of the Vessal potential for types u and

v.

A,, 1. A parameter in the Born-Mayer-Huggins interaction potential between

atoms of type u and v. 2. A parameter in the two-body part of the Vessal

potential for types u and v.

A,, A parameter in the three-body part of the Vessal potential corresponding

to the interaction between types u, v and w.

B 1. A parameter in the Born potential. 2. A parameter in one form of the

Lennard-Jones potential. 3. A parameter in the two-body portion of the

Stillinger-Weber potential. 4. A parameter in a stocastic crack simulation.

B(9Oy An angularly dependent function used in the Vessal potential.

BI., A parameter in the two-body part of the Vessal potential for types u and

V.

B,, The strength of the three-body interaction of types u, v and w in the

Vashishta potential.

C 1. A constant the maximum crack velocity equation. 2. A parameter in

a stocastic crack simulation.

C11, A parameter in the two-body part of the Vessal potential for types u and

v.

C,, 1. A parameter in the Born-Mayer-Huggins interaction potential between

atoms of type u and v. 2. A parameter in the two-body part of the Vessal

potential.

D 1. Fractal dimension. This is an approximation of the box counting

dimension when applied to real materials. 2. A parameter in a stocastic

crack simulation.

D" Fractional part of a fractal dimension.

D,,, A parameter in the two-body part of the Vessal potential for types u and

V.

D,, A parameter in the Bom-Mayer-Huggins interaction potential between

atoms of type u and v.

E,, A parameter in the two-body part of the Vessal potential for types u and

v.

E Young's modulus or elastic modulus.

EK System kinetic energy.

Ep System potential energy.

F A mathematical set of points.

F(r) Force function of r, based on a potential function of r.

F,,, A parameter in the two-body part of the Vessal potential for types u and

V.

F, The interatomic force vector between atoms i and j.

Fqk The interatomic force vector resulting from the interaction between atoms

i, j and k.

G Shear modulus.

H,, The strength of the steric repulsion in the interaction of types u and v in

the Vashishta potential.

Kic Fracture toughness, the stress intensity required to cause type I fracture.

L 1. A linear dimension that is smaller than an edge dimension of the

molecular dynamics cell, used in the truncation function of the Born-

Mayer-Huggins potential based on the Ewald sum. The value of the force

cut-off distance is often used. 2. Longitudinal modulus.

M, A material dependent constant in the equation relating the failure stress of

a material and the radius at which a transition in the fracture surface

occurs, such as mirror to mist, mist to hackle, or mist to crack branching.

The subscriptj corresponds to the particular transition.

N The number of atoms in a simulation.

N, The number of atoms of type i.

N,(F) A function that gives the number of boxes of side length 6 that cover the

set F.

P System pressure.

P_ Pressure tensor in the mn coordinates.

Pi-, A parameter in the two-body part of the Vessal potential for types u and

V.

Q,,, A parameter in the two-body part of the Vessal potential for types u and

v.

R, A random number evenly distributed between 0 and 1 used in assigning

random velocities to atoms.

R,,, A parameter in the two-body part of the Vessal potential for types u and

V.

R? The crystallographic R factor, a figure of merit used in the comparison of

pair correlation functions where better fits are represented by lower

values.

Sl, A parameter in the two-body part of the Vessal potential for types u and

v.

T Temperature of the system.

T(r) The truncation function applied to the force function in the Soules

potential.

Tp(r,) The experimentally derived pair-correlation function at discrete values of

r.

T,,,(r,) The simulation derived pair-correlation function at discrete values of r.

xvii

V Volume of the molecular dynamics cell.

Z, The effective charge of an ion of type u in the Vashishta and Born-Mayer-

Huggins potentials.

a 1. A parameter in the Morse potential. 2. A cut-off distance used in the

Stillinger-Weber potential.

ao 1. The dimension of a unit cell along the crystallographic a direction. 2.

A parameter having dimensions of length in the relation between fracture

toughness and surface fractal dimension of for a material.

aY The effect of a defect at location (i,j) in a stocastic crack simulation.

b A parameter in the Born-Mayer-Huggins potential.

c Crack length.

co The dimension of a unit cell along the crystallographic c direction.

c. An integer indicating how many cell lengths a particle should be translated

in the u dimension to bring it back into the primary cell.

d, The deviation caused by a defect at location (i,j) in a stocastic crack

simulation.

dj The u coordinate difference between atoms i andj.

d' The u coordinate difference between atoms i and j after correction for

periodic boundary conditions.

dimF The box-counting dimension of set F.

f2(r) A dimensionless two-body potential, used in the Stillinger-Weber

potential.

fo(r1,r2,0) A dimensionless angularly dependent three-body potential, used in the

Stillinger-Weber potential.

gdev(1) A function that returns a random Gaussian deviate with a variance of 1.

i 1. An index indicating an atom. 2. An index representing a location on

a discrete x-y plane in a stocastic crack simulation.

j 1. An index indicating an atom. 2. An index representing a location on

a discrete x-y plane in a stocastic crack simulation.

k An index indicating an atom.

k, Boltzmann's constant.

k,, A parameter in the three-body part of the Vessal potential corresponding

to the three-body spring constant for the interaction between types u, v

and w.

I An unspecified constant in the Vashishta potential.

m, The mass of atom i.

m. The mass of an atom of type u.

m An index representing a coordinate in space.

n 1. A parameter in the Born potential. 2. An exponent in the truncation

function used in the Soules force function. 3. An index representing a

coordinate in space.

n, The number of valence electrons for atom type u, used in the Born-

Mayer-Huggins potential.

p A parameter in the two-body portion of the Stillinger-Weber potential.

q A parameter in the two-body portion of the Stillinger-Weber potential.

q. The charge on atoms of type u.

r Length or distance, often corresponding to an interatomic distance.

ro The cutoff distance for the three-body interaction in the Vashishta

potential.

r4, The decay length of the r" interaction term in the Vashishta potential.

r, The position vector for atom i.

ro, The position vector for atom i at time t.

rv The interatomic distance between atoms i andj.

Radius about the initiating flaw where a transition in the fracture surface

occurs, such as mirror to mist, mist to hackle, or mist to crack branching.

The subscriptj corresponds to the particular transition.

The interaction cutoff distance used in the Soules force function.

The u dimension of the MD cell.

1. An index indicating an atom type. 2. An index representing a

coordinate in space.

1. An index indicating an atom type. 2. An index representing a

coordinate in space.

The velocity of atom i.

The maximum crack propagation velocity in a material.

1. An index indicating an atom type. 2. An index representing a

coordinate in space.

The x coordinate location on a discrete x-y plane in a stocastic crack

simulation.

The x coordinate distance between atoms i and j.

The u Cartesian coordinate of atom i.

The y coordinate location on a discrete x-y plane in a stocastic crack

simulation.

The y coordinate distance between atoms i and j.

The z coordinate of a stocastic crack at coordinates (x,y).

The z coordinate distance between atoms i and j.

The length of a time step.

The electronic polarizability of an ion of type u in the Vashishta potential.

1. The Pauling coefficient, used in the Born-Mayer-Huggins potential, for

atom types u and v. 2. A parameter used in the screening function for the

long range interactions in the Feuston-Garofalini potential between atom

types u and v.

6 A small length.

e 1. Strain. In this work we use engineering strain. 2. A parameter used in

Lennard-Jones and Stillinger-Weber potentials to scale the potential

energy. 3. A constant in the truncation term of the Born-Mayer-Huggins

potential function based on the Ewald sum, typically between 0.175 and

0.35.

Y 1. Surface energy. 2. A parameter in the three-body portion of the

Stillinger-Weber potential.

X A parameter specifying the strength of the three-body interaction in the

Stillinger-Weber potential.

X,, A parameter that gives the strength of the three-body interaction for atom

types u, v and w in the Feuston-Garofalini potential.

,1- The exponent of the steric repulsion in the interaction of types u and v in

the Vashishta potential.

o0 A parameter in the Morse potential.

42,,(r) Pair interaction potential function of r for atom types u and v. The

subscript 2 indicates that this is a two-body interaction potential. If only

one atom type is present, the subscripts uv are not written.

43.,-(r,, r,, r Three-body interaction potential function of positions r,, r, and r, for atom

types u, v and w. The subscript 3 indicates that this is a three-body

interaction potential. If only one atom type is present, the subscripts uvw

are not written.

O1.w(ri, rIk,' 0,) Angularly dependent three-body potential of r,, r,, and the

included angle, 0k for atom types u, v and w. The subscript 0

indicates that this is a simple angularly dependent three-body

potential. If only one atom type is present, the subscripts uvw are

not written.

p 1. Density. 2. Hardness parameter in the Born-Mayer-Huggins potential.

a 1. Stress. 2. A parameter used in Lennard-Jones and Stillinger-Weber

potentials to scale the interaction distances.

ac Critical stress required to cause fracture in a sample with surface flaws.

of Failure stress.

aTh Theoretical cohesive strength.

a, Radius repulsive parameter for atom type u in the Bom-Mayer-Huggins

potential.

00 The tetrahedral bond angle, cos-'(-1/3), in the three-body portion of the

Vessal potential.

0,k The i-j-k bond angle.

0,, The preferred bond angle for the three-body interaction of types u, v and

w in the Vashishta potential.

r7, Shear in the u-v plane: displacement in u dependent on change in v.

v Poisson's ratio.

xxii

Abstract of Dissertation Presented to the Graduate School

of the University of Florida in Partial Fulfillment of the

Requirements for the Degree of Doctor of Philosophy

ATOMIC-SCALE DYNAMIC PROCESSES

IN THE BRITTLE FRACTURE OF SILICA

By

THOMAS P. SWILER

April 1994

Chairperson: Joseph H. Simmons

Major Department: Materials Science and Engineering

The molecular dynamics (MD) simulation technique was used to study the atomic-

scale dynamics processes that may take place in the fracture of brittle materials. Silica

was chosen as the host material for this study because it occurs in both crystalline and

vitreous phases, it is of great technical interest, and it has been widely simulated using

MD.

Two models for silica based on different interatomic potential functions were used

and compared. The Soules two-body potential resulted in structures having better

densities and closer correspondence to neutron diffraction data than did the Feuston-

Garofalini three-body potential, which resulted in fewer topological defects and a

cristobalite structure that did not undergo a symmetry transformation under uniaxial

strain. Both vitreous silica and cristobalite were modeled.

Fracture experiments were performed by applying uniaxial strain to the various

samples at strain rates ranging from 0.05 ps' to 5.0 ps'. The fracture behaviors for all

samples were dependent on the strain rate. Samples strained at higher strain rates

exhibited higher fracture strengths, less localized fracture, and rougher fracture surfaces.

Cristobalite samples were stronger than vitreous silica samples.

Both permanent structural variations and thermal vibrations affected the path of

fracture through the samples. The fracture path was extremely sensitive to initial

conditions, indicating that fracture is a truly chaotic process.

We found that the effects of the strain-rate dependence could be traced to the

"flow" of strain throughout the sample during fracture. In a normal fracture process, the

local strain-rate at the crack tip may increase as a crack propagates, causing fracture to

appear as a high strain-rate process in real material. A fracture mechanism is proposed

that may result in a fracture surface containing self-similar features like those observed

in the fracture of real brittle materials. A simple stochastic crack growth simulation is

performed to show the viability of the proposed mechanism.

xxiv

CHAPTER 1

INTRODUCTION

What is Fracture

The word fracture originated from the latin wordfrangere, which simply means,

to break. In days of antiquity, brittle objects that were whole simply broke into parts if

they were not strong enough to support the applied load. Of course, objects still break

if they are not strong enough to support the applied load, but we now have the tools to

investigate how they break. Fracture is not considered an instantaneous event by which

one piece becomes two or more pieces, but a transformation process that can be followed

from start to finish. However, as much as fracture has been experienced, the process

is not fully understood.

We propose that a complete understanding of fracture can be obtained only if we

understand what happens in a material under strain to fracture at the atomic scale.

Fracture displays a complexity that cannot be fully modeled using a continuum model for

a material body. Griffith' and Orowan2 proposed criteria for predicting the growth of

cracks in brittle materials using strain-energy release and stress-intensity criteria,

respectively, in a body that is assumed to be composed of a continuum of material. The

formalism put forth by these models is still dominant in fracture mechanics. The stress

intensity criterion has been extended to predict slow crack growth of materials in a

1

2

corrosive environment by Michalske and Freiman.3 However, although these models

predict the growth of cracks in an object under an applied load and the role of cracks in

reducing the strength of objects, they cannot predict the complex geometries that the

fracture process creates. Thus, we must advance beyond a continuum model and view

fracture as a process that occurs in a many-body system, instead of as a process that

occurs in a continuum of matter, to understand all the processes that result in the details

observed in fracture.

We seek to model fracture on the atomic scale and then attempt to extrapolate the

behaviors that we observe to the mesoscopic scale.

Fracture as a Transformation

Fracture is a simple process in some respects, and it is nearly infinitely complex

in others.

In terms of fracture, materials fall into two gross categories--brittle and ductile.

Brittle materials are those materials that fracture with little deformation before the

initiation of fracture. Ductile materials are those that are capable of undergoing

significant deformation before fracture occurs. Certain materials may be either brittle

or ductile under different conditions.

Within the field of materials science and engineering, the costs associated with

both fracture and corrosion to the national economy on a yearly basis are often cited to

emphasize the negative aspects of these processes. These are the processes that cause

3

the failure of useful objects. As a result, much of the research in the area of fracture has

been in involved in predicting the strength of materials.

The processes by which the usable strengths of materials in real-world

applications are limited are well understood. The strengths of ductile materials are

limited by the formation and movement of dislocations that cause permanent deformation

at the yield stress. This may cause the parts made from these materials to fail in their

purposes, whether the parts have completely fractured or not. The strengths of brittle

materials are usually limited by flaws that result in regions of significantly higher than

average stresses within the materials, thus causing fracture to be initiated even when the

average stress on the material is less than its cohesive strength.

How do ductile materials fracture when dislocations are pinned, and how do

brittle materials fracture when there are no flaws? Dislocations in engineering ductile

materials are routinely pinned by work hardening and by precipitation hardening so that

processes giving ductility are deactivated and failure occurs by brittle fracture. When

there are no flaws present in brittle materials, the fracture of these materials becomes

dependent on the dynamic aspects of materials.

The importance of the dynamics of failure in ductile materials has long been

accepted. As a result, the dynamics of failure in ductile materials has been extensively

studied. However, the dynamics of failure in brittle materials is not so well understood,

because brittleness has traditionally implied that fracture is complete almost as soon as

it has begun.

In order to study the dynamics of fracture in brittle materials, we require high-

speed investigative techniques on the atomic scale. These techniques do not currently

exist. Thus, in order to study the dynamic processes occurring in brittle fracture, we

perform simulations to examine the behaviors that we cannot see by normal experimental

means.

Molecular Dynamics

Molecular dynamics (MD) is a numeric technique that builds on a knowledge of

interatomic potential functions to investigate the behavior of a system of atoms under

various conditions. It was first used by Rahman4 to investigate the motion of Ar atoms

in a liquid, and has since been used in studies ranging from enzymatic reactions to the

fracture of materials.

Molecular dynamics can be compared with another commonly used simulation

technique, Monte Carlo simulations, by examining the path that the system takes while

being simulated. Monte Carlo simulations make random perturbations of a system to

find favorable configurations. It is not necessary that the individual perturbations be

realistic because only the thermodynamics and structure of the final state are important.

In contrast, molecular dynamics simulations allow movement of atoms according to

Newton's equations of motion. Thus, the Monte Carlo simulation technique allows the

simulation of physically realistic structures, but not physically realistic dynamics.

Molecular dynamics, on the other hand, is capable of simulating dynamics, but not

necessarily the most optimal structures. In MD, movement to optimal structures may be

5

blocked by either inertial barriers, where particles may possess too much kinetic energy

to stay in the deepest potential wells, or potential barriers, where particles may not have

enough kinetic energy to get out of shallow potential wells, within the short time allowed

by the simulation.

Although the limitations of molecular dynamics are numerous, its strengths are

obvious. Molecular dynamics is a technique that can be used to study the atomic

processes and the resulting properties of material systems that are not conducive to a

study by analytic means. A quick inspection of the disordered surfaces of fragments of

a fracture event illustrates the complexities involved in fracture. Thus, molecular

dynamics offers the ability to witness the events that occur during fracture that result in

these surface morphologies.

Processes in noncrystalline solids such as vitreous silica cannot be understood

without such a simulation technique because of the complexities involved. Traditionally,

the properties of glasses have been described by thermodynamics because of an inability

to understand the dynamics involved. Most glasses have an enthalpy function of

temperature that displays a glass transition region. The glass transition region is the

temperature region wherein the enthalpy in a super-cooled liquid ceases to behave as a

liquid and begins to behave as a solid upon cooling. The laws of thermodynamics tell

us that within the glass transition region there must be certain limits on the changes in

enthalpy. For example, since both the enthalpy and the specific heat of a liquid are

greater than those for a crystal, a linear extrapolation of the enthalpy function for a liquid

to lower temperatures may cross the enthalpy function for the crystalline phase.

However, because of the excess configurational potential energy of a liquid as compared

to a crystal, the enthalpy of a liquid should never be lower than the enthalpy of a crystal

at the same temperature. This is the Kauzman Paradox. In a sense, it guarantees that

the liquid of the composition must go through a glass transition if it does not crystallize

first. After a material passes through a glass transition upon cooling, its structure

becomes configurationally arrested, and both glass and crystal phases have nearly the

same specific heats. At this point the glass and crystal enthalpy curves become parallel,

so that the glass curve will not cross the crystal curve, therefore avoiding the paradoxical

situation. Simulation is important in the investigation of these transformations because

one can pinpoint the types of liquid-like behaviors possible at different temperatures,

rather than lumping behaviors into a featureless normal distribution. Thus,

nonequilibrium molecular dynamics has become an important technique in the study of

glasses.

It is in this backdrop that simulations of fracture in a brittle material, a glass,

vitreous silica are performed. Molecular dynamics not only gives us the ability to

simulate a dynamics event such as fracture, it provides us also with a way to obtain the

structure, in the case of a glass, so that these simulations can be performed on facsimiles

of real materials. Limited as it may be, MD gives us the ability to study the events that

we cannot view directly.

Statement of Problem

Vitreous silica is a material in which the intrinsic fracture behavior has become

important. Silica optical fibers are routinely made to near theoretical strength by

eliminating surface and bulk flaws. Despite this control over the state of flaws in this

material, variations in the strength of samples still exist. The Griffith criterion, when

applied to the variation in sample strengths, predicts flaw sizes on the order of a few

angstroms. This work seeks to find the source of this variation in strength in this

material.

Additionally, the dynamics of fracture, although not directly observable in fracture

experiments, leave a signature on the resulting fracture surfaces. Reading this signature

is the subject of the science of fractography. Fractography has been successful in

locating the flaw that initiated fracture and in finding the state of stress in the material

when fracture was initiated. However, the dynamic processes that result in the exact

path that fracture takes after being initiated have remained elusive. The simulation

technique that we use should be able to shed some light on these processes.

CHAPTER 2

BACKGROUND

In this work, we use molecular dynamics simulations to study the atomic-scale

fracture processes in silica. Therefore, we will review the important aspects of this

study. Since we hope to generalize the results of our brittle fracture simulations to brittle

fracture in general, we will start with a discussion of fracture. Then we will examine

the various phases of silica that are at the center of our particular study. We will then

discuss the technique of molecular dynamics, including the limitations of this technique

as well as the extent of applications. Various simulation models for silica will then be

reviewed. We will then describe how fracture is simulated and studied using molecular

dynamics. Finally, we will identify the aspects of fracture that a dynamical technique

such as MD may allow us to study. Because of the various fields that are drawn

together, a smooth progression between topics will not be attempted in the background

section.

Brittle Fracture

Characteristics of Brittle Fracture

Brittle fracture is the mechanical failure of a material without plastic deformation.

Brittle fracture is not limited to traditionally brittle materials. Metals or polymer

materials that are considered to be plastic or ductile can be brittle under certain

conditions of temperature and strain rate. Even inorganic glasses that are the epitome

of brittle solids at room temperatures behave as viscous liquids at temperatures above the

glass transition and are not brittle at low strain rates.

Experimentally, brittle fracture is generally a sudden process, so the dynamics of

this process are not well known. Observations made during the fracture process by

Dickinson et al.5 show that brittle fracture is accompanied by an ejection of particles and

light. To date, measurements of these fractoemitted particles remain one of the few

methods for characterizing the dynamic processes leading to brittle fracture.

The Ultimate Strength of Materials

The ultimate strength of a material without pre-existing surface flaws that fails by

brittle fracture can be calculated based on a knowledge of its structure and the properties

of its interatomic bonds. This value is known as the cohesive strength of the material,

and it can be calculated for a material given the Young's modulus, the surface energies,

the atomic structure, and an assumption of the form of the interatomic potential.267

10

A cohesive strength calculation may use one of many pair potential functions.

Two possible potentials are the Morse potential and the Born potential. The Morse

potential is applicable to covalent systems and has the form 4(r) =4e 'r-r) -2e -(r-

where 42(r) is the pair potential energy function of r, r is the interatomic distance, ro is

the equilibrium interatomic distance, and &, and a are constants that can be found from

materials properties. The Born potential is applicable to ionic systems and has the form

02(r) =A +B where 4 and r are as described above, and A, B, and n are constants that

r r"

can be found from materials properties.

Cohesive strengths are derived from these pair potential functions by calculating

the derivatives of the functions and finding an expression for the maximum stress

encountered in breaking the chemical bonds. Stress is defined as .= 1 d (r) where a

A dr

is the stress, and A is the cross-sectional area covered by the bond. The maximum stress

is then found where i(r) is a maximum, or when d 2(r) is zero. Thus, the

dr dr2

theoretical cohesive strength is found by calculating the value of r for which d ,2(r)

dr2

is zero, then substituting this value of r into the expression for 1 d0(r)

A dr

The resulting expressions for the cohesive strength for the two potentials are, for

the Morse potential, % Yy and, for the Born potential, E[ 2 r) where y

is the surface energy, E is the Young's modulus, and n=-rE-

For vitreous silica, with E = 7.3 x 01O Pa, y = 4.56 J/m2, and r. = 1.62 x 10-10

m, the theoretical strength given by the Morse potential is 23 GPa. This compares

favorably with the highest strength of 14.3 GPa measured for treated silica fibers in

liquid nitrogen by Proctor et al.,' but is about two orders of magnitude greater than the

measured strength for bulk silica or untreated silica fibers.

The Effect of Flaws

Discrepancies between theoretical strengths and measured strengths are common

in brittle materials. This occurs because the presence of flaws in materials result in local

strains that are much greater than the average strain within the material under stress.

Plastic deformation does not readily occur in brittle materials, so the stresses in these

areas of locally high strain cannot be relieved. These regions of high stress cause failure

to occur in a material even when the average stress throughout is much less than the

theoretical cohesive strength. Because of this, flaws may be considered as stress

concentrators in the fracture of brittle materials.

12

Griffith' calculated the tensile stress, a,, which could be supported by a sample

of a material containing a crack of a given length, c. His calculation was based on the

assumption that a crack will propagate under an applied stress when the rate of decrease

in elastic strain energy during crack propagation is greater than the rate of increase in

surface energy resulting from the formation of new surfaces. Thus, for a sheet with an

edge crack of length c, the critical stress is given by' -_ 2Ey Similar expressions

Irc

have been derived for various types of crack geometries.

Regimes of Brittle Fracture

The limiting processes in the brittle fracture of glasses in reactive environments

are different in the three regimes categorized by Wiederhorn.' Fracture processes of

regimes I and II are limited by environmental effects and occur at stresses, or stress

intensities, well below the theoretical strength of the material, and crack growth

velocities are affected by species in the chemical environment. In type III fracture, the

material becomes unstable under the applied stress without regard to environmental

effects and the high crack growth velocities that result are limited only by the mechanical

properties of the material.

Environmental Effects

Stress-induced corrosion or slow crack growth corresponding to types I and II

fracture are important in limiting the strengths of materials in many applications.

1 r--. T -T

In water

0.1 In moist nitrogen

In dry nitrogen

0.01

0.001

0.0001

I e-05

le-0 -

le-09

0.2 0.3 0.4 0.5 0.6 0.7 0.6

Stress Intensity, MPa'm-0.5

Figure 2-1. Schematic showing crack propagation velocities as a function of stress

intensity factor. Adapted from Michalske."

Michalske" studied the fracture surfaces created by slow crack growth to better

understand the mechanisms involved. In this type of fracture, environmental species

diffuse to the crack tip, weakening the chemical bonds at the tip and causing crack

growth at stress intensities less than o, but greater than some minimum stress. The

velocity of the crack tip is then controlled by different mechanisms in the two regions.

Since the reactivity of the material to corrosion at the crack tip is dependent on the local

stress, the rate of reaction and therefore the crack growth velocity is dependent on the

14

stress intensity at the crack tip in type I fracture. The rate of diffusion of the

environmental species to the crack tip controls the crack growth velocity and is nearly

independent of the stress intensity in type II fracture. Thus, the crack growth velocity

in region II has a lesser dependence on stress intensity than in region 1. A schematic

showing this behavior is given in Figure 2-1.

Crack Prooaeation Velocities

The maximum velocity of crack fronts in type III fracture is limited by the speed

at which information of the approaching crack front may travel through the system. In

a mechanical system, this speed is the speed of sound. Different methods have been

applied to calculate to maximum crack propagation speed in real materials, giving similar

expressions of the form2 v=C E where C is a constant between 0.4 and 0.6, and p

Fm\ P

is the density of the material. Maximum crack velocities are found by experiment to be

consistently less than, but on the same order as, the values given by the above

expression.

Reading the Fracture Surface

Because brittle fracture is a sudden process that does not easily lend itself to real-

time study, much of the analysis of brittle fracture is performed by fractography, the

study of fracture surfaces. Since we seek to describe the processes occurring during

fracture, we need to understand what fracture surfaces can tell us about these processes.

Through fractography, one can determine the location of the flaw that initiated fracture,

determine the size of the initiating flaw, and determine the stress at which the sample

fractured by examining the patterns on the fracture surface.

Brittle fracture surfaces have four distinct regions--the smooth mirror region near

the initiating flaw, the rougher mist region beyond the mirror region followed by the

coarse hackle region, and finally the crack branching region extending beyond the hackle

to the limits of the fractured sample. By locating the smoother regions of a fracture

surface, one approaches the area of the initiating flaw, and by advancing to the center

of the smoothest region, the mirror region, one can pinpoint this flaw.

Fracture surfaces can provide information on the strain state of the material before

fracture as well as the location of the initiating flaw. The transitions between the mirror

and mist regions, the mist and hackle regions, and the hackle and crack-branching

regions are distinct. It has been found that the radii of these regions are related to the

failure stress by the relationship13 oaJ=M. where a, is the failure stress, r. is the

transition radius corresponding to either the mirror-mist, the mist-hackle, or the hackle-

crack-branching transitions, and M, is a material dependent constant for each of the

transitions.

Recently efforts have been made to describe the aspects of the fracture surface

that have traditionally eluded a meaningful quantitative description. These aspects are

the seemingly random fluctuations that exist within all regions of fracture on varying

16

scales. Fractal geometry has been used to describe the statistical nature of these

fluctuations, as will be discussed later in this review.

Silica

Bonding in Silica

Silica is the chemical compound SiO2. Oxygen is a group VI element and silicon

is a group IV element. Oxygen is more electronegative than is silicon, but the difference

in electronegativities is not great enough to make the bonding in SiO, purely ionic, so

the bonding has significant covalent nature. Because silicon is a group IV element, it

requires four bonds to complete its valence if covalent bonding is assumed, or it has four

valence electrons to donate if ionic bonding is assumed. Oxygen is a group VI element,

so it requires two bonds to complete its valence if covalent bonding is assumed, or it will

gain two electrons to complete its valence if ionic bonding is assumed.

Structure of Silica

The basic structural unit of most condensed silica phases is the silica tetrahedron

as shown in Figure 2-2.14 In each silica tetrahedron, there is one silicon atom

tetrahedrally coordinated by four oxygen atoms. Each oxygen atom is then shared by

two such tetrahedra to maintain local charge neutrality. Because the silica tetrahedron

is the basic structural unit of silica glass as well as most crystalline phases of silica, even

the glass phase is not without short-range order on the scale of a silica tetrahedron.

Because of this, we refer to this phase as vitreous rather than amorphous.

(a)

(b )

Figure 2-2. The silica tetrahedron is the basic structural unit of silica structures.

a) One silica tetrahedron. b) Two tetrahedra linked by corer sharing.

0 is the Si-O-Si bond angle. Adapted from Holloway.'4

Phases of Silica

Because silica structures are relatively open, they are relatively free to change

symmetry with temperature. Thus, silica comes in several polymorphic phases. The

major crystalline phases of silica are quartz, the phase thermodynamically stable at

temperatures below 1140 K, tridymite, stable in the temperature range of 1140 K to 1770

K, cristobalite, stable between 1740 K and the melting point of 2000 K, and coesite, the

high pressure phase.'5.'6 In addition, there are further distinctions within the individual

phases. The quartz, tridymite, and cristobalite phases are subdivided into high and low

subphases. The transformations between the subphases are displacive, meaning that there

is no rearrangement of the structural elements necessary for transformation, only

reorientation of the structural units. The transformations between the major phases are

reconstructive, meaning that there are structural rearrangements required to go from one

phase to another. The vitreous phase is the usual result of cooling molten silica to a

solid when no effort is made to induce crystallization. For a summary of the phases of

silica, see the review article by T.D. Coyle."

In crystalline silica, the tetrahedral units are linked in an orderly fashion, resulting

in a structure possessing long-range order. In vitreous silica, the tetrahedra are linked

in an apparent disordered fashion. This results in a structure possessing short-range

order given by the rigid constraints on the size and shapes of the tetrahedral units, which

diminishes rapidly beyond the first coordination distance because of the less ordered

Table 2-1

Characteristics of the Phases of Silica

Phase Symmetry Density (g/cm3) Tet Linkage

Low Quartz Hexagonal 2.6 Corner

Tridymite Hexagonal 2.22 Corner

Cristobalite Tetragonal-Cubic 2.33 Corer

Coesite Monoclinic 3.01 Corner

Keatite Tetragonal Corner

Shishovite Tetragonal 4.35 Octahedrally

Coordinated

Silica-W Tetragonal 1.97 Edge

Vitreous NA 2.20 Corner

linkage of the tetrahedral units. The characteristics of the phases of silica are given in

Table 2-1.

The Crystal Structure of Cristobalite

Cristobalite is of interest to us because its high symmetry, low density structure

is probably the most comparable to the vitreous phase. At room temperatures, both

phases are essentially thermally arrested high temperature phases. Thus cristobalite is

the crystal phase of choice to compare to vitreous silica in studies of the effects of

structural order and disorder in silica.

Although quartz may be the equilibrium phase, the cristobalite phase can be

relatively stable when quenched to low temperatures because the reconstructive

20

transformations between the major phases is slow. Within the cristobalite polymorph,

two sub-polymorphs are present, high, or beta, and low, or alpha, cristobalite, which

undergo a displacive transformation within the temperature range of 470 540 K. This

displacive transformation results in a symmetry change from the high temperature cubic

symmetry to the lower temperature tetragonal symmetry upon cooling.

The structure of high cristobalite"1 is similar to the diamond structure. The silicon

atoms are in the diamond positions and oxygen atoms are located between each pair of

silicon nearest neighbors. The symmetry of high cristobalite is cubic, and the structure

belongs to the Fd3m (Oh7) space group. The unit cell of high cristobalite has an edge

length of 7.16 A at 560 K and contains eight formula units.

The structure of low cristobalite has the same connectivity as the high cristobalite

structure. Although this phase has tetragonal symmetry with prism dimensions

ao = 4.96 A and co = 6.92 A, the structure is pseudocubic when the cell is rotated by

450 about the c-axis and the sides a. are transformed into diagonals of the base of the

prism, that is, ao = 2a, = 7.02A and co = co 6.92A. In this description, the

structure is a relatively small distortion of the high cristobalite structure. The unit cell

of low cristobalite contains four formula units and belongs to the space group P412,

(D,4).

The Rine Structure of Silica

A useful description of intermediate structure in the phases of silica is the ring-

size distribution. The statistical nature of this measure allows it to be used for vitreous

silica structures as well as crystalline structures. The ring structure plays prominently in

the description of silicate structures because these structures are relatively open as a

result of the low coordination of silicon and oxygen atoms, allowing identification of

individual rings in the network.

A ring in a silica structure is a closed loop of adjoining tetrahedra that are

identified by the number of tetrahedra contained in the ring rather than by the number

of atoms that are bonded to form the ring. Work by Galeener" has shown that different

numbers of tetrahedra in rings give rise to differences in the energies of the structure.

Works by Marians and Hobbs20,21 and Marians and Burdett22 have described the

topological constraints on the ring size distributions in crystalline and vitreous forms of

silica and on vitreous silica surfaces. Most crystalline forms of silica are composed

entirely of six-membered rings, whereas vitreous forms necessarily have some ring sizes

different from six.

The ability to make vitreous silica by a sol-gel route has enabled a detailed study

of the role of rings in the formation of a silica network structure. In sol-gel processing

of silica, silicon alkoxide precursors in an aqueous solution undergo condensation

reactions such as the following, where alkoxides are identified by the formula unit

(OR):23

Si(OR)4 + H20 (RO)3Si-OH + ROH

2 (RO),Si-OH 4 (RO)3Si-O-Si(OR)3 + H20

(RO)3 Si-OH + Si(OR)4 7 (RO)3Si-O-Si(OR)3 + ROH

22

As the polymerization process continues, a network structure gradually forms. Because

of the thermal and environmental constraints on the system, the resulting gel is not

generally structurally similar to bulk silica. The absence of a high temperature step in

the gelation process prevents the system from exploring different structural configurations

that would result in the system seeking a global metastable state similar to that obtained

by fused silica. Instead, the system structure is the result of a growth process by which

Si(OR)4 molecules react first to form chains, then to form a network structure as

SiO,(OR)(4 groups become linked with an associated loss of water molecules. Because

the interim structure results from the dynamic growth process emanating from reaction

centers, the structure tends to be fractallike--the geometry of the structure grows from

a center point and through a series of branches, attempts to fill space. Silica gels are

transformed into bulk silica structures by the processes of drying, stabilizing, and

densifying. During drying, liquid water and reaction by-products are removed from the

gel. In stabilization, physically adsorbed water is removed. Finally in densification,

chemically adsorbed water is removed, and the structure attains full connectivity. During

the densification step, structural changes occur that elucidate some of the medium-range

structural details of silica.

Raman studies of gels during the densification process tell an interesting story.

Superimposed on a somewhat broad Raman spectrum are narrow "defect" bands at 490

and 608 cm', so called because their presence cannot be explained by continuous random

network models.24 As gels are stabilized by driving off chemically combined water in

the structures and then densified by collapsing the porous gel structure, the intensity of

23

these defect bands increase, then decrease. This has permitted identification of these

bands with the condensation of dangling silanol groups on the surfaces of the gels. As

these dangling bonds become linked, highly strained three- and four- membered rings are

formed and subsequently destroyed as the gel densifies and the internal surfaces are

eliminated. This shows that the concept of rings corresponds to a real structural

component in silica glass.

Elastic Properties of Vitreous Silica

Mallinder and Proctor" found that vitreous silica at low temperatures displayed

a nonlinear elastic modulus in which the apparent elastic modulus is increased at high

strains. At 77 K, the Young's modulus, E, and the shear modulus, G, as functions of

strain, E, were found to be E=71.9(1 +5.75e) GPa and G=31.5(1+3.06e) GPa,

respectively. The authors attributed this behavior to the opening of Si-O-Si bond angles

in the structure, causing the elongation of silica rings. They reasoned that extension of

the structure could be more easily accomplished by initially opening the Si-O-Si bond

angle than by elongation of Si-O bonds.

The same authors also found that sodium silicate glasses exhibited a nonlinear

elastic modulus in which the apparent elastic modulus decreased with increasing strain

under similar conditions. In this case the authors reasoned that the extra cations were

positioned in the structure so as to prevent the rotation of tetrahedra under low strain.

The function for the elastic modulus given for this material was given as

E=71.1(1-5.11e) GPa.

24

Thomas and Brown26 studied the elastic properties of fused silica and found that

this material exhibited a similar nonlinear elastic modulus at room temperature. The

authors suggested that this anomaly may play an important role in crack propagation and

fracture.

Proctor2 examined the effects of a nonlinear elastic modulus on the fracture

behavior of fused silica. The author concluded that the nonlinear elastic effect must be

considered when calculating either failure strains in high strength samples or stress values

from failure strains, but that the nonlinear effects on the calculations of the size of

fundamental Griffith flaws are self-cancelling and generally can be neglected.

Kulawansa et al.2 performed an atomic force microscopy (AFM) study on silica

fracture surfaces. The authors of this study analyzed the details of the fracture surface.

They explained that some of the surface features observed were the result of a crack

blunting mechanism caused by a hardening of the material under high strain, as they

claimed would be consistent with a nonlinear elastic modulus that increases with strain.

The Fracture Strength of Silica

Measurements of the strengths of silica fibers indicate that nearly flaw-free fibers

have been produced. Proctor et al.8 performed strength measurements on as-drawn fibers

and found a narrow distribution of material strengths, indicating that samples were nearly

flaw free. When the samples were fractured in liquid nitrogen, the authors obtained a

maximum strength of 14.3 GPa. Kurkjian and Paek9 performed an analysis of the

tensile strength statistics for silica fibers and found that most of the variability of the

25

apparent strengths of fibers could be attributed to the variability of fiber diameters along

the gauge length of the tensile test, giving a nearly single-valued strength of 5.75 GPa

with a variation of less than 1%.

In a seminar at the University of Florida Department of Materials Science and

Engineering in the Spring of 1990, Kurkjian noted that measurements suggested that fiber

strengths were not completely single valued and indicated the presence of flaw sizes on

the order of a few angstroms. A scanning tunnelling microscope (STM) was used to

examine the topography of as-drawn fiber surfaces to explore the nature of the possible

flaws.

Proctor27 interpreted fracture strengths for silica fibers and found that average

stresses were well below the highest measured value. Using elastic modulus data

obtained in earlier work, he obtained values of flaw sizes of 3.6 nm for fused silica

fibers and 7 nm for soda-lime glass. He then made the assumption that the most

prevalent surface flaw is the largest broken ring of silica tetrahedra on the surface of the

fiber and thus deduced average ring sizes of 5 to 7 tetrahedra for fused silica and 8 to

12 tetrahedra for borosilicate glass. The corrected strengths obtained from the later

strength data ranged from 17.9 to 18.7 GPa for fused silica and 6.8 to 10.1 GPa for

borosilicate glass.

Stress Corrosion of Silica

A model for the chemical reactions resulting in the stress corrosion of vitreous

silica was proposed by Michalske and Freiman3 and was further explored by Michalske

0

0 o

so d H^

S HH

% O O SO

HHH

0o oW 0 0 0

(a) (c) (c)

Figure 2-3. Chemical reaction proposed by Michalske and Freiman for

environmentally enhanced Si-O-Si bond rupture. The strained crack tip

bond (a) adsorbs a water molecule, (b) is cleaved by dissociative reaction,

and (c) is converted to surface silanol groups. Adapted from Michalske

and Bunker.30

and Bunker.30 In this model, the authors propose that a water molecule aligns itself with

a strained Si-O-Si bond, resulting in a weakening of the bond, and the water molecule

is consumed in a reaction that breaks the bond and creates a pair of silanol groups. By

this model, three properties may be common to all active corrosive species: (1) the

molecule has at least one lone (nonbonding) electron pair (Lewis base), (2) the molecule

possesses a labile proton (Bronsted acid), and (3) the distance between the acid and base

sites conforms with the Si-O bond distance (0.163 nm). A schematic of this process is

given in Figure 2-3 and shows three stages in bond cleavage.

Hench and West''32 studied the reaction paths in the fracture of Si-O bonds in

rings of various sizes using molecular orbital calculations. They found that when

particular bonds of unconstrained rings were strained, the energies required to break the

bonds were significantly lower than the expected energy of the Si-O bond. This behavior

resulted because unconstrained rings of various sizes, each having one strained Si-O

bond, underwent a ring contraction mechanism whereby another tetrahedron within the

ring approached the oxygen atom of the strained bond, five-coordinating the tetrahedron,

and weakening the initially strained Si-O bond. The authors found that 4-membered

rings were the weakest rings, with an activation energy for Si-O bond breakage of 320

kJ-mol', followed by 3-membered rings with an activation energy of 400 kJmol', and

5-membered rings with an activation energy of 430 kJ-mol"'. When water was

introduced into the system, the various sized rings showed different tendencies for

hydrolysis, the critical step in environmentally assisted fracture. In this case, the authors

found that 3-membered rings were the easiest to hydrolyze, with an activation energy of

30 kJImol'l, followed by 4-membered rings, with an activation energy of 120 kJImol'

and 5-membered rings with an activation energy of 160 kJ-mol-. In both the water-free

and water-present cases, the authors concluded that the ring structure of bulk silica

influences the fracture of Si-O bonds, and must be considered when studying this

process.

Because water is an active component in the stress corrosion of silica, and since

it is prevalent throughout the atmosphere, it plays an important role in all fracture

experiments on silica except for those experiments in which it is purposely and carefully

excluded. Understanding the role of water is of vital importance when interpreting the

observations of any fracture experiment. The fracture of silica in the absence of water

is quite different from the fracture of silica in a water-free environment, as demonstrated

by a comparison of strengths measured by Proctor8 (14.1 GPa) in a water-free

environment at 77 K and by Kurkjian and Paek29 (5.75 GPa) in a laboratory environment

at room temperature.

All simulations to be presented in this work pertain solely to a water-free

environment.

Molecular Dynamics

Survey of Applications

A molecular dynamic (MD) simulation is, as its name implies, a technique that

models the movement (the dynamics) of particles (of molecules). As pointed out by

Hoover3 in an introductory article on the subject, interesting results can be obtained from

even the simplest of systems because of the many-body interactions that may be involved.

Molecular dynamic simulations are not difficult to visualize. Imagine the

collective motion of several balls on a billiard table with a frictionless surface and

perfectly elastic bumpers. Imagine now a computer model of this motion, and you have

envisioned an MD simulation of a two-dimensional hard-sphere system with rather

peculiar boundary conditions. Further, MD simulations are capable of extracting

thermodynamic quantities such as pressure and temperature from the dynamic system,

29

providing a simulation link between the system of discrete particles in motion and the

thermodynamic quantities that describe that system.

Molecular dynamic simulations were first used by Rahman4 to study the

correlations of atomic motions in liquid argon. MD was first used in simulations of

liquids for various reasons. Clearly, liquids are dominated by particle dynamics,

therefore, their properties can be quite well simulated by MD methods. Since liquids are

fluid, their structures are sufficiently dynamic so that slow nucleation and growth

processes are not required to obtain reasonable simulations of their structures, and

comparatively short MD simulations are capable of simulating reasonable structures.

Finally, despite the dynamic nature of liquid, there are strong interactions among atoms

and transient structures of sufficient interest to warrant a simulation of structure and

properties, as opposed to an ideal gas where there is little interaction between atoms and

no structure.

The applications for MD simulations have expanded over time. Non-equilibrium

MD is concerned with the dynamics of systems and not in thermodynamic equilibrium,

and is used in simulating glass structures. The structure and dynamics of materials

surfaces have been studied. For example, Heyes et al.3435 and Heyes3" studied the

structure and dynamics of surface and near-surface regions in ionic crystal films.

Broughton and Gilmer37n."39 studied crystal-fluid interfaces of Lennard-Jones metals.

Gilmer, Grabow, and Bakker4 used MD simulations to model epitaxial growth in

Lennard-Jones systems. Briinger et al.41 reported the use of MD simulations in the

30

determination of crystallographic structures of protein molecules from X-ray and neutron

diffraction data.

There has been recent interest in using MD to study the interactions between

surfaces in a field called "nanotribology." Landman and Luedtke42 reported MD

simulation results in studies of the interaction of metal tips with metal substrates.

Harrison et al.43 reported the use of MD in studies of friction between diamond surfaces.

Belak et al." used MD simulations to model atomic scale deformation of surfaces by a

process similar to machining.

The evolution of computer hardware, along with the evolution of computer codes,

is continually increasing the capabilities of MD simulations. MD simulation codes are

capable of being run efficiently on parallel computers with some modification.

Furthermore, researchers have developed specialized computers for performing certain

simulations more quickly and efficiently. An example of such a computer is described

by Bakker et al.45 A result of the development in computer hardware is that systems

composed of tens of millions of particles can now be simulated on parallel computers

consisting of thousands of processors, whereas systems having fewer than one thousand

particles were simulated in the infancy of MD.

Recently, MD simulations have reached popular appeal with the advent of

commercially available MD simulation software packages, such as the Insight package

from BioSym Inc. One can now buy a simulation system that is capable of simulating

a wide range of materials and configurations. It is likely that the use of MD in all types

of materials simulations will increase dramatically since it is no longer necessary to either

write simulation codes or modify existing simulation codes for a wide range of

applications. Despite the advent of commercially available simulation codes, many

researchers will most likely continue using custom-made codes, which are optimized for

the simulation of certain systems or to run on specific computer hardware.

Survey of Potential Functions

The identity of a simulation is strongly dependent on the potential function used,

which in turn is based on the type of bonding in the material being simulated. There are

several common potential functions that are used in MD simulations.

In order to simplify the description of the various potential functions in common

usage, we will use the same means of summation for all potentials, and common symbols

whenever possible to permit comparison between the potentials (see the Key to Symbols

in the Appendix). Summations are not necessarily listed as they would be performed the

most efficiently in the computer code, but are written for logical simplicity. The general

interaction potential for a system will be written as

where N is the number of atoms in the system,

i, j, and k represent individual atoms as specified in the summations,

u, v, and w represent the atom types of i, j, and k, respectively,

2l.(r0,) is the function representing the two body portion of the potential for types

u and v, and

32

-3,_(r,, r,, r) is the function representing the three body portion of the potential

for types u, v, and w. Since the three body portion of the potential is

most often a function of two bond lengths and a bond angle, it may

usually be equivalently referred to as .ik(rj, rj, OJ). We have written

the summations to work most directly with this form because it is

applicable among all the three-body potential functions that we will

review.

The use of indices u, v, and w to denote atom types is a departure from the usual

notation in the literature. This notation is used to unambiguously differentiate the

constants or functions are dependent on atom types only from the variables that are

dependent on individual atoms. Although this is usually clear from context, there have

been some cases of ambiguities and errors in the literature. For systems of only one

atom type, the subscripts u, v, and w will not be used.

A common pair-potential used for the condensed phases of noble elements is the

Lennard-Jones 6-12 potential. The function for this potential is A2(r)=--_ A+ where

r r12

A and B are positive constants and r is the distance between atoms. This potential is

usually written in the dimensionally more appealing form (r)=4ef2(), where

f2(r) = r12 r' is a dimensionless function, and a = (BIA)"6 and e = A2/4B are

constants.46 This potential was also commonly used in the simulation of metals, since

it often results in a close-packed crystal structure similar to the structure of many metals.

33

Recently, a potential based on the embedded-atom-method (EAM) developed by

Daw and Baskes47 has become a popular method of simulating processes in metals with

empty or filled d bands.4" This potential uses the assumption that atoms are embedded

into a local electron density provided by the remaining atoms of the system. The

electron density function of the system is a many-body interaction problem, so this

potential is a many-body potential. This potential has an advantage in metals because it

is less empirical than is the Lennard-Jones 6-12 potential traditionally used in these

systems.

Covalently bonded materials require a potential with an angularly dependent term

to simulate the angularly dependent bonds that result from bonding through molecular

orbitals. One such potential was developed for silicon by Stillinger and Weber.49 An

angular dependence is obtained by adding an appropriate three-body term to a truncated

Lennard-Jones-type potential. The form of the potential is 02(r) =f' for the two-

body part, and (r, r 9~ )e.r' ( k o) for the three-body part. The function f in the

two-body part of the potential is

[A(Br P-r )exp( r

(r) r-a

0, r>a

while the function f, in the three-body part of the potential is

Xexp (Y Y coso+ I2,
f0(r,r6k j rq
[0, r >a V r >a,

and the constants e and a serve to make2 and f dimensionless functions.

Ionic materials are generally modeled using a potential that incorporates a

coulombic potential with an exponential repulsive term. The Born-Mayer-Huggins

potential is such a potential. It has the general form

(.) +A=- e brc rj S+D rS

where r> is the inter-ionic distance, q, and q, are the ionic charges, and A,,, b, C,, and

D,, are constants.

Potentials for Silica

There are several commonly used potentials for modeling silica. The first of

these is a variant of the Born-Mayer-Huggins described by Busing,so first used in silica

by Woodcock et al." and later modified by Soules. The second is a variant of the Born-

Mayer-Huggins potential that incorporates a three-body term similar to that used in the

Stillinger-Weber potential. These two potentials will be described in detail in the

procedure section of Chapter 4.

Vashishta et al.52 derived a potential function that included a three-body part and

applied it to various silica polymorphs. The two-body portion of the potential is given

by

H qq (aZ,2 Z2)

S(r r. 2i

and the three-body portion of the potential is given by

B ,exp (-- -- coos ,- r)
0,._(r,rO = k) r-ro i -roj

0, r> ro V > r

with constants H,, and 7,, representing the strengths and exponents of the steric

repulsion, Z, and a. representing the effective charge and the electronic polarizability

of ions of type u, r4, representing the decay length of the r' interaction term, B,,, ro,

and 0,, representing the strength, the cutoff distance and the preferred bond angle for

the three-body interaction, and I an unspecified constant. The authors claim that this

potential gives the correct Si-O-Si bond angle distribution where earlier potentials failed.

Another potential function that is now commonly used for silica is described by

Vessal et al.53 This potential function has the most adjustable parameters of any

reviewed here, and claims an excellent fit to experimentally measured properties of

various phases of silica. Wright54 found this potential created a vitreous silica structure

36

having the best fit to neutron diffraction data of any MD simulated structure at the time

of the comparison. The authors also claim agreement with the presumed Si-O-Si bond

angle distribution in real silica. The functional form of this potential is quite complex.

The two-body portion of the potential has long-range and short-range parts. The long-

range Coulombic portion of the potential is calculated using the Ewald55 method applied

to a fully ionic system of silicon and oxygen. The short-range part of the potential is

given by the four-part Buckingham potential

Aexp(r', r

A ..r +B l.u +Cl.r3 +Dr42? +El.r, +Fl.., r,,,< r < r2.

P.u pr+Qiuri2 +Rui +Sr, r2~

C,

-- r3., < < r.

where r, is the interatomic distance between i andj, and most other constants are purely

empirical. The three-body potential for O-Si-O interactions has the form

U(rrk,Ok) = iA (B(GO) )'exp(- r )

where A, is a constant given by A k B(O,) is a function given by

2(00-r)2

B(O,) =(00-7)2-(,, -I)2, and the constants kk and 00 are the three-body spring constant

and the tetrahedral bond angle, respectively. With various spline approximations, etc.,

this potential has more than 40 parameters. With so many adjustable parameters, it is

not surprising that this potential function can claim the best fit to the experimental

properties of silica, although the simultaneous optimization of all parameters is

impressive nonetheless.

Alavi et al." proposed a potential that adjusts the effective charge of silicon and

oxygen ions based on their coordination. This potential was proposed to correct

problems created when network-modifying ions are included in the structure in which

some of the network forming ions are not fully coordinated. The concept of variable

charge may be important in fracture studies as well where coordination of atoms change

during fracture. Tsai et al." have used such a mechanism in the simulation of fracture

in silicon. Thus, this charge-transfer potential may be one to watch in the future.

Applications of MD in Models of Silica

Woodcock et al." first simulated a silica-like structure using a variant of the

Bom-Mayer-Huggins potential as used for ionic systems. These authors acknowledged

38

that a purely ionic potential was not fully applicable to a system with a fair degree of

covalent character, but they tried the potential since it was available and easy to apply.

Despite the known limitations of this potential, comparisons of simulation derived pair

correlation functions (PCF's) with PCF's deduced from x-ray diffraction data showed

rough agreement and a comparison of system energies showed that the potential energy

calculated from the simulation was 5% less than an experimentally derived Born-Haber

lattice energy for quartz. The authors made observations and performed experiments on

the simulated system. They measured diffusion constants at various temperatures and

pressures, and studied the permanent compaction of silica by the application of a

compressive stress. Mitra et al.58 performed simulations of vitreous silica structures

using a variant of the Bom-Mayer-Huggins potential and reported good agreement to x-

ray and neutron diffraction data. These authors also examined the internal energy as a

function of temperature and reported that the simulated system exhibited a glass transition

at 1500 K.

Soules59 performed a simulation of sodium silicate glasses using a potential similar

to Woodcock et al.5" Again, the structure was reported to be in significant agreement

to x-ray diffraction data. Soules found that alkali and alkaline earth cations clustered

around non-bridging oxygens, as is generally predicted for these structures. Coordination

information showed that the number of non-bridging oxygen atoms was strongly

dependent on the amount of alkali ion present. Ring size distribution data showed

changes with varying alkali content as well, although it was not clear from the limited

extent of the reported ring size distribution what the changes were. Soules and Busbey6

39

further performed simulations to study diffusion of sodium in alkali silicate glasses and

found that diffusion coefficients obtained from simulation at high temperatures were in

agreement with extrapolated experimentally measured values. They then concluded that

purely ionic potentials for silica were sufficient to model transport properties.

Garofalini and Levine," and Feuston and Garofalini62 examined the structural

changes on surfaces of vitreous silica. These simulations created arbitrary surfaces by

removing periodic boundary conditions between a pair of sides of the MD cell. The

general finding of these simulations was that less-highly charged ions tended to segregate

toward the surfaces. In pure vitreous silica, these species were oxygen ions, while in the

alkali silicates, alkali ions also tended to segregate toward the surfaces.

Garofalini" studied the reactivity of water on vitreous silica surfaces. He found

that water reacts with a silica surface through a series of penta-coordinated silicon

transition geometries in which water molecules are somewhat paired on the surfaces.

Comparisons of Simulated Structure to Experiment

Researchers generally do not like to point out the discrepancies of their models

in the published results of their MD simulations. They often compare radial distribution

functions derived from simulation with those deduced from diffraction experiments, and

if there are the same number of peaks in the distribution function and they occur in the

same general location, they call it a "good" fit. However, based on some minimal

constraints on the structure such as the preferred Si-O bond length for silica, in no case

can radial distributions be too far wrong. Wright"" has attacked the casual comparisons

40

of RDF's derived from MD simulated structures to the RDF's deduced from

experimental techniques such as neutron diffraction data and has proposed a standard for

comparing the experimental and simulated RDF's. This standard is a variant of the

crystallographic R factor, and is given by R = Tp( Ti (r) By this standard,

5 Ttp(r,)

the current champion in the fit of simulated structure to experiment is Vessal et al." with

an R7 factor of 9.0% .

MD Simulations of Fracture

The use of molecular dynamics for studying the failure of materials is not new.

Two approaches have been taken in these studies. The first has been to perform fracture

experiments in simplified MD systems in order to perform detailed studies of the

elements of fracture. The second has been to study fracture using MD simulations that

attempt to model real systems.

Paskin et al.,6,"7 Dienes and Paskin," and Sieradzki et al.6 performed simulations

of crack growth in two-dimensional crystalline lattices. These works demonstrated the

applicability of the Griffith criteria in simple 2-D systems modeled by molecular

dynamics, and studied mechanisms of crack growth, including the formation of

dislocations and speed of crack propagation, as functions of applied stress, crack length,

41

and potential parameters. These works provide an atomistic basis to fracture mechanics,

however, they do not provide details on the fracture of particular materials.

MD simulations of fracture in silicate glasses were first performed by Soules and

Busbey70 on sodium silicate glass samples at elevated temperatures. The periodic

boundary conditions used made this experiment equivalent to the fracture of a fiber with

a diameter of about 50 A, so this study did not model fracture in a bulk material. The

essential characteristic of glass fibers in the real world that gives them interesting

mechanical properties is the absence of surface flaws, resulting in nearly theoretical

strengths, but the simulated "fibers" did not share this characteristic. At the atomic

scale, the surfaces appear rough and would not be considered pristine if the simulated

"fibers" were scaled to the size of real fibers. Thus, behavior similar to fibers is not to

be expected, and is not observed. What is observed is that samples undergo necking as

the simulated surfaces dominate their mechanical properties. Cavitation within the

samples is also observed, indicating that the system sometimes takes this route to absorb

the added volume while maintaining connectivity and minimizing the total energy at each

stage of strain.

Soules71 examined the role of the vitreous structure in the fracture of silica by

comparing fracture in vitreous silica to fracture in high cristobalite silica. He

investigated how the Si-O-Si bond must break during fracture and found that for the bond

to break under an applied strain, the bridging oxygen must vibrate out of its potential

well. Thus, he showed that fracture in silica is an activated process and undergoes a

relaxation. He also found that high cristobalite exhibited a high strength of 70 GPa as

compared to a strength of 24 GPa for vitreous silica.

Kieffer and Angell72 applied an isotropic tensile strain to a silica sample in order

to rupture the structure. They found that the resulting structures had a fractal character

similar to that of aerogels, and that the fractal dimension was dependent on density but

independent of strain rate, indicating that the fractal dimension for silica under a given

set of thermodynamic constraints is unique.

Ochoa and Simmons7 examined the role of strain rate in the fracture behavior of

vitreous silica. They found that within the range of 0.4 4.0 ps', the strength, as

defined by the maximum stress supported by the system during strain to failure,

increased as a function of strain rate from about 45 GPa to 95 GPa. They attributed the

decrease in strength with decreasing strain rate to a relaxation process where atomic

vibrations caused some reorientation of the structure capable of lowering its cohesive

strength. They also found that the effects of strain on the structure were dependent on

strain rate, with the highest strain rates causing an elongation of Si-O bonds and lower

strain rates resulting in an opening of Si-O-Si bond angles in the strained structures prior

to fracture. Finally, they discovered that the maximum stress did not occur at a strain

immediately prior to rupture of the structure, but rather at some point where some

fraction of the bonds were broken internal to the structure.

Ochoa et al.74 examined the role of thermal vibrations in the fracture of vitreous

silica and high cristobalite by simulating the fracture of these materials at strain rates

ranging from 0.2 8.0 ps-'. They used a potential function based on the Born-Mayer-

43

Huggins potential with two sets of potential parameters to investigate the sensitivity of

results on a specific potential. The authors found that the maximum stress during strain-

to-fracture increased with increasing strain rate for the vitreous silica samples, ranging

from 30 GPa to 65 GPa, but was nearly independent of strain for the cristobalite sample

at 80 GPa. They explained that the decrease in strength of the amorphous structure with

decreasing strain rate was a result of rotation of (SiO)4' tetrahedra during strain-to-

fracture at low strain rates that was not allowed at high strain rates. They found that

voids placed in the samples resulted in a decrease in the strengths of the samples and

altered the mechanism of fracture at low strain rates. Finally, they found that the

differences in the potential parameters did not significantly affect the behaviors observed.

Simmons et al.75 examined the mechanisms required for stabilization of silica

surfaces following fracture. They found that fracture is stabilized by the movement of

dangling silicon atoms back into the structure, leaving predominantly oxygens on the

surface, resulting in a repulsion between the two freshly created fracture surfaces.

Chaos in Fracture: Determinism and Fractal Surfaces

We earlier discussed brittle fracture in silica. That discussion was based on the

principles of fracture mechanics--maximum strengths of static structures, the effects of

flaws on the strength of materials, and the role of environmental species during fracture.

However, although fracture mechanics may serve to identify the causes of fracture, it is

ultimately unable to predict the path of fracture that results in the seemingly random

structure of fracture surfaces. The complexities of fracture surfaces indicate that the

44

particular dynamics of the fracture process is important in determining the path of

fracture. Although molecular dynamics simulations may be capable of elucidating the

role of dynamic processes during fracture, we cannot analyze every atomic motion in the

fracture of a simulated system. Because we can easily observe the deviation of the

system when we make a small perturbation in the starting conditions, the concept of

chaos with its sensitivity on initial conditions gives us a handle on analyzing the complex

behavior that cannot be analyzed on an atom-by-atom basis. Thus, without looking at

all dynamic interactions in detail, we can find the significance of these interactions in a

qualitative way.

In nature, fractal-like shapes often result from chaotic processes. This is seen in

brittle fracture as well. Thus we will also discuss fractal geometry and its appearance

in brittle fracture.

Investigations of Chaotic Processes

"Chaotic" is a descriptive term applied to systems that follow paths which are

strongly dependent on initial conditions. The term, "strongly dependent on initial

conditions" implies that any deviation in initial conditions ultimately results in system

paths as divergent as system paths resulting from large deviations in initial conditions.

The concept of chaos originated from simulation work of weather patterns by

Lorenz.76,77 Lorenz found that small deviations in the starting conditions of a numerically

solved differential equation resulted in a divergence of the system behavior over time.

Although the paths of the systems remained close for some time, they gradually deviated

to the point that they no longer possessed any correlation that would indicate a shared

beginning. The paths taken by the simulated systems deviated as a result of the loss of

precision associated with converting the computer's internal representation of numbers

into a decimal representation of limited precision, then reconverting this representation

into an internal binary representation when re-starting the simulation. Although the

simulation was not necessarily capable of a fully accurate representation of the weather,

it demonstrated that without a total knowledge of starting conditions and external effects,

the prediction of weather for indefinite periods of time would be impossible. The general

implication of Lorenz's work is that the behavior of systems described by coupled

nonlinear equations cannot be predicted within a certain range in the long term. In an

MD simulation where a material described by coupled nonlinear potential functions,

similar behaviors may be expected.

Mathematical descriptions of chaos are particular. One of the tools used in

investigating the chaotic nature of a system is a Poincard map, which is a map of the path

of a system in a space of some or all of its variables. Chaotic systems are often found

to have maps that approach but never achieve a periodic behavior, and such a map was

constructed by Lorenz. Often a projection of these maps creates either a fractal or an

intricate pattern. Another method used in investigating the chaotic nature of a system

is an investigation of an apparently randomly varying quantity. In fractoemission

experiments, Langford et al.78 performed autocorrelation, Fourier transform, and box-

counting dimension analyses on the curves of fractoemission intensity as functions of

time. They found evidence of deterministic chaos from these analyses, and related it to

the fractal dimension of the resulting fracture surfaces as found by the slit-island

technique. We will discuss the implications of a fractional dimension in a real material

later in this review.

Predictability of Brittle Fracture

Fracture occurs when a mechanical system becomes unstable by the application

of strain. What does it mean when a system becomes unstable? A system is unstable

when it is in a state that will undergo a transformation to a favored state under any

perturbation. This means that as a system becomes unstable, even previously

insignificant dynamics within the system may play a role in determining the path to the

favored state. This seems to be the case in brittle fracture where the fracture mechanics

that can quite accurately predict when a system will become unstable is ultimately unable

to predict the precise path to fracture after it has become unstable.

The choices available in predicting the path of brittle fracture are quite limited.

Either brittle fracture either occurs by cleavage, where fracture surfaces approach atomic

smoothness, or it does not. In cleavage, the path to fracture is known or controlled. But

in all other fracture events, the path is not known and is not exactly controlled. For

example, distributions of sample strengths have not been completely explained by known

sample flaws as pointed out by Kurkjian previously. Also, according to a private

discussion with T.A. Michalske in 1991, there appears to be some uncertainty in the

dynamics of slow crack growth in pristine silica under vacuum where there should be no

environmental effects, and thus, no effect of stress on the rate of crack growth. The

47

question arises--is the inability to completely describe the fracture process the result of

this process being chaotic?

Fractal Geometry

A fractal is a mathematical set of points that is best described by fractal geometry

and has a fractional dimension. Fractal geometry is a branch of mathematics that

encompasses Euclidean geometry as a special case--the case of integral dimensions.

Objects in geometry are idealized and are not realized in the physical world, and are

referred to as sets. Although a set is a mathematical abstraction, physical objects may

approximate sets. Thus, when a surface is described as being fractal, it really means that

the shape of the surface approximates a fractal set, just as when a physical object is

described as spherical, the description implies that the shape of the object approximates

a sphere to certain tolerances. Thus, physical objects may befractal-like, but they can

not be fractal.

One of the most fundamental descriptions of a set is its size. The measurement

of sets is made using mathematically rigorous techniques. In order to measure a set, the

proper dimension for the set must be used. The Hausdorff dimension for a set is the

dimension in which a meaningful measurement of the set can be made. As its definition

is beyond the scope of this discussion, refer to Falconer79 for a complete definition of the

Hausdorff dimension.

The Hausdorff dimension is the only dimension for which the set may have a

finite non-zero measure. For example, in Euclidean geometry, the measure of an area

48

is infinite if the measurement is made using a measure of length, is zero if the

measurement is made using the measurement of volume, and can only be non-zero and

finite only if the measurement is made using the same dimension of the set, that is, using

a measurement of area. However, in contrast to examples of sets in Euclidean geometry,

it is not always obvious in what dimension a measurement of a non-Euclidean set should

be made. Thus, although the Hausdorff dimension is not in itself a measure, it is

essential to the measurement process, and it is a useful descriptor for a set.

The Hausdorff dimension is the most fundamental definition of the dimension of

a fractal set, but other definitions are also used. The box counting dimension operates

on the graph of a set, and is defined as dim F lim logN(F) where F is the set that

-o -log

describes the object, the function N, gives the number of boxes of linear dimension 8

required to cover F.79 In many cases, the box-counting dimension is equivalent to the

Hausdorff dimension.

Fractal sets are of interest to mathematicians because of their unusual properties,

but why are sets that are best described using fractal geometry of interest to physical

scientists? The answer is that shapes found in nature often approximate the shapes of

fractal sets. Since mathematicians have discovered the generating functions for creating

fractal sets, physical scientists may assume that the physical processes that result in the

formation of such shapes are described by similar functions.

Several authors have related how physical objects approximate fractal shapes.

The "proof" that they give to show that these objects are fractal is that the "fractal

dimension," a limited approximation to the box-counting dimension, is not an integer.

The definition of "fractal dimension" used by Mecholsky et al.80 is given by the relation

Nr = 1, where D is the "fractal dimension," N is the number of boxes of size

proportional to r that cover a set, or the number of ruler lengths proportional to r that

measure the length of a curve. These are related to quantities in the definition of the

box-counting dimension where dimB F -D, N, N, and 5 r. An important distinction

is that no limit is involved in the definition of the "fractal dimension" as described here.

Thus, the quantity D is a useful description of an object that approximates a fractal over

a limited range of r. However, it does not prove that the object is fractal, only that it

possess some properties that approximate those of a fractal object over a limited range

of scale.

Although the Hausdorff dimension remains the most precise definition of the

dimension of a fractal set, the term "fractal dimension" is often used in the physical

sciences to denote the approximate box-counting dimension as defined above. Because

of this, we will follow the same usage for the remainder of this work. However, one

should bear in mind that this fractal dimension is obtained only over a limited range of

scale, and the validity of extrapolations to different size scales must be tested.

Applications of Fractal Geometry to Fracture

Recently, there has been a great deal of interest in describing fracture surfaces

using fractal geometry. There are several reasons for doing this.

50

The first is that all statistically significant details of these seemingly random

surfaces can be described with a limited amount of information, because the same surface

features may be found over a wide range of scales, with the atomic scale being the

absolute lowest limit. Thus, the description of fracture surfaces using fractal geometry

provides a possible link between those observations made by fractography and those

made by atomic scale simulation.

The second is that a scalar quantity, the fractal dimension, is obtained for the

description of fracture surfaces. Just as the pair correlation function is a useful

description of a random glass structure, the fractal dimension is a description of the

random structure of a fracture surface. However, a fractal dimension, like a pair

correlation function, does not unambiguously describe the structure. It is merely a

statistic of the structure that may be compared to test the possible similarities of various

structures. Despite its limitations, the fractal dimension can be related to material

fracture properties, as will be shown.

Finally, and probably most important from the standpoint of understanding

fracture, a description of a fracture surface as fractal-like may help in the understanding

of fracture processes. Fractal sets can be created by iterated generator functions. By

describing a fracture surface as fractal-like, one may be able to find the generating

function for the surface, and then elucidate some of the processes which occur during the

progression of fracture. Since fractal-like shapes often imply a chaotic generation

process, do the fractal characteristics of fracture surfaces imply that fracture is chaotic?

Later in this work, we will attempt to show how fracture may proceed to form fractal-

like surfaces.

The Fractal-like Nature of Fracture Surfaces

If fracture occurred only by the orderly cleavage of crystal planes in crystals and

it always created minimum energy surfaces in glasses consistent with the applied

mechanical constraints, then the science of fracture would be nearly complete. However,

brittle fracture is not that simple. During the formation of fracture surfaces, the effects

of the dynamic lattice seem to counter the attempts of the material to minimize its surface

energy. Even in the fracture of a highly ordered material such as single-crystal silicon,

Tsai and Mecholsky" found that the resulting fracture surfaces were not smooth, but had

a fractal-like character. Mecholsky et al.0O reported that fracture surfaces of ceramic

materials also exhibit fractal-like characteristics from investigations of alumina and glass

ceramic fracture surfaces. Such surfaces indicate that a dynamic process is at work in

determining the path that fracture follows in brittle materials.

Because the description of fracture surfaces as being fractal-like implies that

fracture surfaces may have detail to nearly the atomic scale, such a description may

provide a link between macroscopic measurements and microscopic or atomic scale

behavior. Thus, the finding that fracture surfaces have a fractal-like character opens new

avenues for study into brittle fracture.

Mecholsky et al.82 reported a relationship between the toughness of a material and

the measured fractal dimension of a fracture surface of that material. They found that

52

the toughness was related to the fractional part of the fractal dimension by the relation

Kc =E/oD where Kc is the fracture toughness, E is the Young's modulus, D" is the

fractional part of the fractal dimension, and ao is a constant with units of length, which

may be representative of some characteristic length in the material or of the process zone

during fracture. Since surfaces with a higher fractal dimension generally appear more

tortuous, an easy way to remember this relationship is the phrase, "the rougher, the

tougher."

Applicability of Fractal Analyses to MD Simulated Systems

An important question that needs to be answered is the capability of modeling

fractal-like structures by MD simulations. Quality fractal analysis requires a range of

scale that is not readily available to MD simulations. Can one derive any valid results

from an analysis of fracture surfaces formed in MD simulations.

Kieffer and Angell7 reported on the isotropic fracture of silica. Under these

conditions, a bulk silica sample appeared to attain a structure reminiscent of silica

aerogels. By using a box-counting technique applied to pair-correlation functions, they

found the fractal dimension of such samples to be less than two.

In meetings with John J. Mecholsky and Joseph H. Simmons, we discussed what

would be necessary to perform an analysis of fractal dimension on surfaces created by

uniaxial fracture. Mecholsky used a slit-island technique that is essentially a box-

counting applied to a surface on real fracture surfaces. We looked at Kieffer's paper as

53

a model to see how a fractal analysis might be accomplished on the fracture surface

created by an MD simulation. Kieffer's method viewed the fractured silica as a mass

fractal, whereas the method used by Mecholsky examined surface fractal characteristics.

In a mass fractal analysis, an object may have a dimension ranging from three,

corresponding to a continuous three-dimensional network solid, through two,

corresponding to a layer structure, to one, corresponding to a linear molecule. All

samples are composed of a finite number of atoms that can be arranged in three-, two-,

or one-dimensional arrays. Since scattering or pair correlation functions treat atoms as

point particles, this analysis has no bias towards finding a set to be three dimensional,

and the dimension found is wholly dependent on the arrangement of particles in three-,

two-, or one-dimensional structures.

Clearly, this analysis was not suited for finding the surface fractal dimension of

a solid as 2 + D' where 0 < D' < 1. Such an analysis would be best performed by

some kind of slit-island analysis. There are two problems associated this type of

technique. The first is that it is not a simple matter to define the fracture surface in a

MD system because the system is logically a system of point particles. The second is

that our samples are too small to perform a fractal analysis, because samples of size on

the order of tens of Angstroms and composed of units on the order of Angstroms could

not have the range of scale required to perform an accurate analysis of fractal dimension.

Thus, although a slit-island analysis could be performed in an MD system, it might not

give an accurate representation of the true fractal character of the surface.

Y.L. Tsai et al." recently conducted a molecular dynamics study of the fracture

of single crystal silicon on the atomic scale. This system was chosen because of previous

measurements of the fractal dimension on fracture surfaces of this material using a

scanning tunneling microscope. Also, this system had already been simulated using

molecular dynamics, and the potential function (Stillinger-Weber) was similar to one that

was already written into code (Feuston-Garofalini). The authors were successful in

performing a slit-island analysis of a fracture surface simulated by MD, and were able

to calculate a fractal dimension, although the range of measurement on which this

dimension was based was necessarily limited.

The technique used by Tsai differs from the technique used by Kieffer in an

important respect. Tsai's method is dependent on the space occupied by the constituent

atoms of a material, while Kieffer's method is dependent only on the spatial distribution

of atom centers. Tsai assumes that the constituent atoms can be described as spheres of

a convenient radius while Kieffer assumes that the constituent atoms are point particles.

This difference in approach is due to the difference between surface and mass fractals.

An important effect that results from this difference is that the method that used by Tsai

is dependent on the shape of the space chosen to envelope the atom center, or simply,

the shape of the atom. Since atoms do not have unambiguous shapes, the possible effects

of the shape chosen on the dimension measured is worrisome. On the other hand, the

creation of a surface by carefully chosen atom shapes may compare well to the surface

that may be measured by a scanning tunneling microscope or an atomic force microscope

if such techniques were truly capable of single atom resolution. A fractal dimension

55

obtained by this technique would then be comparable to one obtained from a microscopic

analysis on real fracture surfaces. Results of that study show a strong correlation

between fractal dimensions and surface roughness calculated from the MD simulations

and those calculated from experiments for 100, 110, and 111 surfaces.

Because of the problems associated with finding the fractal dimension of a

fracture surface created by MD simulations of fracture, we will not report any such

analysis in this work. We will, however, try to relate the atomic processes we observe

to possible generating functions that may result in the formation of fracture surfaces

having a fractal-like character.

CHAPTER 3

PROCEDURE

Usine Molecular Dynamics to Model Fracture

In this work, we conduct a series of simulations of brittle fracture in a simple

"model" material for which a wealth of experimental measurements is available.

Although limitations of the modeled interatomic potentials and the spatial and temporal

extent of the simulated systems make it impossible to supply an exact match of the

experimental strength of the material or to calculate exact crack growth velocities, the

simulations were performed to investigate the underlying mechanisms for the trends

observed in fracture.

One of the biggest challenges in modeling fracture on the atomic scale in real

materials, and especially in glasses, is sorting out all the interacting processes that occur

during fracture. Logically, a material may be treated as a large system of interacting

particles. The dynamics of such a system are difficult to predict, and the cause and

effect relationships are difficult to describe analytically. Often the symmetry of such

systems can simplify the analysis, but this approach can not been applied to problems of

fracture where the order of systems is broken. In some cases, the modeling of fracture

in symmetric systems can be treated analytically, as in modeling cleavage fracture where

the system symmetry is dominant over the system dynamics. When a system has little

57

structural order at the outset, as in the case of a glass, an analytic treatment fails, and

a numerical treatment must be used.

Authors who present results from MD simulations often try to directly relate their

results to experimentally obtained results. An example of this is the direct comparison

of the system pressure as a function strain in MD systems to stress-strain curves of real

materials undergoing fracture. However, molecular dynamics gives a perspective of

materials fracture on a microscopic scale that cannot be directly compared to materials

fracture on a macroscopic scale. Therefore, the results of MD simulations are best used

to build an understanding of the processes that occur on the microscopic scale. With an

understanding of microscopic fracture processes and how such microscopic processes

might effect fracture on the macroscopic, one may be able to predict what the

macroscopic experimentally observed behavior will be, possibly to explain previously

unexplainable behaviors. It is not reasonable to expect that MD simulations will give

results which will directly relate to the macroscopic.

The potential function gives the simulated material its molecular identity.

Boundary conditions, system size, system density (for constant volume simulations), and

temperature serve to describe the equilibrium state of a system. In non-equilibrium

molecular dynamics, the path by which the structure is derived is also significant, and

needs to be studied. The study of fracture using MD then breaks down into a simulation

of the sample structure followed by the simulation of strain application and examination

of the resulting processes.

Potential Functions

Modeling of physical systems is possible because many of the interactions within

a physical system can be ignored and fairly accurate predictions can still be made about

certain aspects of the system. Molecular dynamics modeling is possible in our systems,

for example, because the collective behavior of a large group of atoms is relatively

independent of the motion of electrons about individual atoms. In an MD system, the

processes that one intends to model are atomic motions that are controlled by interatomic

potential functions. Because MD approximates electronic bonding structure with

interatomic potentials and ignores electronic states, it cannot be applied to study the

formation of electronic defects in materials, but it has shown great success in studies of

multiple particle dynamics.

Although the fracture of a bond begins with a restructuring of the electronic

configurations forming the bond, the role of nearest neighbors becomes more significant

as the bond is stretched and, as will be shown, dominates the bond fracture process. The

complications caused by environmental species on bonding are eliminated by examining

fracture in a vacuum. In this case, again as will be shown later, MD calculations yield

a bond fracture process similar to the molecular orbital models proposed by West and

Hench." The MD simulations further benefit from the ability to examine properly the

effect of the broken bond on the behavior of the ensemble of neighboring atoms. A

problem not examined here is that of environmentally driven fracture, because of the

major role played by the chemical interaction of environmental species (H' and OH in

water) on the strained bond. In this latter case, ab-initio approaches that examine the

59

details of the electronic configurations of the bond are essential, as shown by West and

Hench."

It is recognized that the forces between atoms are the result of the interaction of

electrons about the individual atoms, however, good results can be obtained by assuming

that the electrons are in some sort of stationary configuration, and that interatomic forces

can be modeled from some pseudo-potential. Systems of atoms with spherically

symmetric electron configurations are well modeled by a two-body potential, which

represents purely ionic character. Systems of atoms with non-spherically symmetric

electron configurations are better modeled by a three-body potential. A three body

potential that models a covalent-like electron configuration thus attempts to model a

covalent system.

In the case of silica, the potential functions that are commonly used either model

silica as a purely ionic system (Woodcock), a pseudo-ionic system (Soules), an ionic

system with a pseudo-covalent character (Feuston-Garofalini), or a pseudo ionic-covalent

system (Vessal). All these potential functions have limitations.

For our simulations of fracture, we have picked two potential functions that have

demonstrated good fits to various pieces of experimental data for silica, and yet are quite

simple. These are the truncated Bom-Mayer-Huggins potential developed by Soules83

and the empirical three-body potential for vitreous silica developed by Feuston and

Garofalini4 that includes an angularly dependent term. The Soules potential is purely

pseudo-ionic in nature, while the Feuston-Garofalini potential incorporates a pseudo-

covalent nature.

The Born-Mayer-Huggins Potential

The basis of both functions is the Born-Mayer-Huggins potential. This potential

function has the form

r) P. b.exp ) + Cr-' + Dr-',

r p

where u, v signify atom types,

0,, is the potential between atoms of type u and v as a function of r,

r is the separation distance between atomic centers,

q,, q, are charges on atoms of types u and v respectively,

e e, are radius repulsive parameters,

b is a positive constant,

p is a hardness parameter,

C,,, D. are the coefficients of dipole induced dipole and dipole -

induced quadrupole terms, respectively,

.,, is the Pauling coefficient, given by

Z Z

0 = (1 + -U + -),

n. n

where Z,, Z, are the units of charge of type u and v, respectively, and

n., n, are the numbers of valence electrons of u and v, respectively.

The multiple terms become important when the atoms are highly polarizable and are

separated by small distances, but are generally neglected in the simulation of silica glass.

Some have questioned the accuracy of this simplification since oxygen exhibits a high

degree of polarization and the atoms are closely spaced in condensed phases. The

justification for ignoring these terms is that atoms within a solid are generally constrained

to the point that these terms would have minimal effect, and these terms involve many-

body interactions that are computationally difficult. When the multiple terms are

neglected, the potential function simplifies to

42(r) = .exp(-- --)+ .

By collecting terms, this becomes

2, (r)A exp=(-) qiq

where

A_, -= b expi .

P

From here, the Soules and the Feuston-Garofalini potential functions take divergent

forms.

Truncation or Screening Functions

In a simulation, it is desirable to have a potential of limited range in order to

minimize the number of particle interactions. However, the Coulombic portion of the

Born-Mayer-Huggins potential approaches zero slowly as the reciprocal of distance and

many neighbors would be required for calculations of all potentially significant

interactions. To overcome this problem, a screening or truncation function is used to

62

limit these interactions. This truncation process simulates the result of a lattice of

alternating charges as used in the calculation of a Madelung constant.

One method of dealing with the slowly convergent Coulomb potential is to

perform an Ewald sum." If one knows the symmetry of a crystal, the electrostatic

potential of ions in the crystal can be more quickly summed by performing the sum in

Fourier space. This is the basis of the Ewald sum approach. From the Ewald sum, a

screening factor is derived that breaks up the Coulombic part of the potential into a local

part that is summed over individual ions, and thus contributes to inter-particle forces, and

a background part that is just added to the overall potential energy and has no effect on

inter-particle forces. Strictly speaking, the results of an Ewald summation are only valid

for the crystal system for which they were derived, however, the functional form of the

screening factor is useful even when the inverse lattice sum has not been performed for

a particular system. Thus, the Ewald convergence parameter may be used in the

summation of the system potential energy using the modified Bom-Mayer-Huggins

potential function

2(r) A_ ,exp(] r q. erfc (,

where erfc() is the error function complement function,

e is a constant, typically 0.175 0.35,6 and

L is a linear dimension of the molecular dynamics cell or r,.

Table 3-1

Parameters Used in the Soules Force and Potential Functions

Potential Function Parameters

p = 0.290023 A

Asisi = 3300 p ergs

As,-0 = 3000 p ergs

AO-O = 2300 p ergs

qs, = +4e

qo = -2e

where e is the electron charge

Truncation Function Parameters

r. = 5.5 A

n=3

The term eL is termed the convergence parameter. It is dependent on an arbitrarily

chosen constant and the size of the molecular dynamics cell, but it is independent of the

type of interaction.

The Soules Potential Function

The Soules force function has been shown to simulate vitreous silica structures

with approximately the correct density and with nearest neighbor distances in significant

agreement with neutron diffraction data. It is the product of the force function derived

from the Born-Mayer-Huggins potential and an empirical truncation function. A direct

differentiation of the Born-Mayer-Huggins potential gives

Fd-, (r) A~ exp )

dr p p r2

200

O-Si -

0-0 ---

Si-Si

50 v

-50 k

-100 1 1 1L

0 1 2 3 4 5 6

Distance, Angstroms

Figure 3-1. Plot of the Soules potential function for various interactions.

The truncation factor has the form T(r) = [ (r )]. The Soules force function is

then

Fo...v(r) = F sv(r) T(r)

S exp(-) + -]. -(

65

The potential function that corresponds to the Soules truncated Born-Mayer-

Huggins force function is

'Sol = 2 Fs^'oe2uv(r) dr

--A J exp(- 1.,]I-(r

i- exp(-'+ t Q -(rr-[ dr,

which is not analytically integrable. Thus, the potential energy of the system is

calculated using a numeric integration of this expression. The parameters used in the

Soules force function are given in Table 3-1. A plot of the Soules potential function is

given in figure 3-1.

The Feuston-Garofalini Potential Function

The three-body potential developed by Feuston and Garofalini is a sum of the

modified Born-Mayer-Huggins potential an empirical convergence parameter and the

three-body portion of the Stillinger-Weber potential.4 The two-body portion of the

potential is

2.(r) = A. exp(- + erfc( )

Notice that the eL term of the convergence parameter has been replaced by /,. This

makes the screening function independent of cell size, but dependent on the type of

200- -- -

O-Si -

0-0 ---

', ', ~ Si-OSio' -- .. ..

150 ;i-Si

100

S50

CL

0 -

a-o

-100

0 1 2 3 4 5 6

Distance, Angstroms

Figure 3-2. Plot of the two-body portion of the Feuston-Garofalini potential function

for various interactions.

interaction, and adds another set of adjustable parameters. The three body portion of the

potential is

h exp[ + osOkCcsO < rf,. A k< r

46.U(r,,k(r -) -r2-) Q .) -rr<)

0, r,_ r V Tjk -> rt

The parameters used in this potential function are given in Table 3-2 and a plot of the

two-body portion of this potential is given in figure 3-2. This shows a shallower

Table 3-2

Parameters used in the Feuston-Garofalini Potential Function

Two-body Potential Parameters

p = 0.290023 A

As,.s, = 1880 p ergs

As,,o = 2960 p ergs

AO-0 = 720 p ergs

qsi = +4e

qo = -2e

where e is the electron charge

3s,i = 2.29

si-o = 2.6

o-o = 2.34

Nonzero Three-body Potential Parameters

Xo-s-o = 180 p ergs

Xsi-o-si = 3 p ergs

mo-si-o = 2.6 A

Ys,-o-s, = 2.0 A

rso = 3.0 A

rsi-o.sC = 2.6 A

cos 0'o-si-o = -1/3

cos 's,.o-s, = -1/3

potential well in the Si-O interaction and a weaker repulsion in the Si-Si and 0-0

interactions than is obtained for the Soules potential.

Simulation Techniques

Integration Algorithms

Molecular dynamics is a computationally intensive technique in which the

simultaneous integration of the motion of many individual particles is performed over

short time steps.

Integration of functions may be performed in many different ways, but only a

couple of different techniques are commonly used when integrating the equations of

motion of a molecular dynamics system because they are well known to provide

reasonably good accuracy with a minimum amount of computation. These are the Verlet

"leapfrog" algorithm and the Gear Predictor-Corrector algorithm.

Verlet algorithm

Our simulations generally made use the Verlet algorithm,5 given by the vector

equation

r,,,. 2. r,, 2 At + FK ,

mi i (j k) i

where r,,,, is the position of the particle after the time step,

r,, is the position of the particle before the time step,

r,,,, is the position of the particle before the last time step,

t is a particular time,

At is the length of a time step,

m is the mass of the particle,

F5 is the sum of all two-body force components on the particle.

F is the sum of all three-body force components on the particle,

where applicable.

Fifth-order Gear predictor-corrector algorithm

The fifth-order Gear predictor-corrector algorithm, as encoded by J.M. Haile,8

was also tested. This was done to verify that the Verlet algorithm, known for its speed

and stability, but not for its overall precision, gives reasonable results as compared to the

slower but more precise Gear algorithm.

The Time Step

The time step is the basic unit of time in a molecular dynamics simulation. The

time step is an approximation for an infinitesimal amount of time that is valid for the

processes under observation. During this fundamental time step, in the Verlet algorithm,

the acceleration is assumed constant. In the Gear predictor-corrector method,

approximations are developed that move the constant function to a higher order time

derivative. In the selection of the fundamental time-step, it is desirable to use a value

short enough to maintain a constant system energy, yet a value long enough to observe

a physically significant process. The usual test for deciding if a time step is of

sufficiently short length is to perform a simulation over about a thousand steps and check

to see if the system energy remains fairly constant. If not, then the time step is too long.

A time step of 0.5 X 10 s was found to be of acceptable length, and was used in these

simulations. This corresponds to about 1/40 of an atomic vibration period.

Limiting Particle Interactions

The magnitude of the task of calculating particle interactions is reduced by

limiting the spatial extent over which these interactions are considered. Neighbors lists

for all particles that include all neighbors which may possibly interact with a central

particle are compiled. Only these neighbors are then considered when calculating

interactions with the central particle. The neighbors lists are updated periodically by

finding all particles within a distance given by the maximum interaction distance plus the

maximum distance a particle may travel during the period between updates. The period

between updates in these simulations is ten time-steps. The maximum distance, found

by Ochoa,8 that a particle may move in a solid in ten time steps at 0.5 X 10"5 s each is

0.5 A. Thus, the maximum distance for inclusion in a neighbors list from the central

atom is 0.5 A greater than the cutoff distance for interaction between particles, or 6.0

A in the case of a interaction cutoff at 5.5 A.

The Pythagorean theorem is used to find the squares of inter-particle distances.

Finding inter-particle distances requires square root operations that are computationally

time consuming, and are avoided. The force table is tabulated in terms of the square of

inter-particle distances, so inter-particle distances are required only when pair correlation

functions or bond angle distribution functions are accumulated.

S* 0* 0 0* 0

r .. m0 ; .

0

0000 0 *0 * *0 00

SPr imary cel

Repea t ed second a r y c e lls

Ne ghborho od abou t t om

Figure 3-3. The effect of periodic boundary conditions in a simulated system. Note

that the neighborhood of an atom may extend beyond the primary cell.

Application of Periodic Boundary Conditions

The use of periodic boundary conditions allows one to model a bulk-like material

with a limited number of particles. A schematic showing the effect of periodic boundary

conditions is presented in figure 2-3.

Periodic boundary conditions are applied to the simulated system when finding

inter-particle distances. If the difference in any coordinate is greater than one-half the

corresponding dimension of the cell, then the corresponding image of the particle located

in an adjacent cell will be closer to the central particle than the particle in the primary

cell. The difference in the coordinate is then adjusted to reflect this by either adding or

subtracting the corresponding dimension of the cell to minimize the absolute values of

the difference. It is the resulting difference that is used in the Pythagorean relation to

find the square of the distance between the particles. The algorithm by which this is

performed is as follows: The differences in coordinates are calculated from

d,i= x, x- ,

and are subsequently adjusted with

S= d., s. round(),

and used in the Pythagorean relation

2 3

where u = 1, 2, 3 denote the Cartesian dimensions of space,

x,, is the u Cartesian coordinate of atom i,

x,5 is the u Cartesian coordinate of atom j,

di is the u coordinate difference between atoms i and j,

d',, is the u coordinate difference between atoms i andj after correction

for periodic boundary conditions,

s, is the u dimension of the cell,

round(x) is a function that returns the integer nearest to its argument, x,

and

r, is the distance between atoms i andj.

The application of periodic boundary conditions gets more difficult when the cell

is sheared. Suppose that one puts the cell oriented with the coordinate axes into simple

shear as shown in figure 2-3. The opposite sides of the cell are shifted with respect to

each other. This must be taken into account when applying periodic boundary conditions

across the shifted faces. The periodic boundary condition equations are

d

c = round( -)

d'l d,, ,, c,) 7.(sC)

d" = d s. round(

u=C.

S-TAl y Sy

---- --------- 0 1----*------ -* -,

/ / / /

*. e-.. ;*. . "

I M - -/

S 0 / *$ e 0

0 ** ** /

Pr ary ce ll

Logical vin ew of slheared primary cell

0

Zanlin Pri may celons

Logical v ew of sheared secondary cells

L___i

Figure 3-4. Schematic showing the application of periodic boundary conditions to a

system under shear strain.

By this method, full periodic boundary conditions are maintained even when the virtual

shape of the cell changes.

Handling Particle Interactions

Interactions between particles are handled by the following process--interaction

geometries are found, interaction forces are found, components of interaction forces are

75

summed on all participating particles, and particle positions are updated to reflect the

effect of interactions.

In the case of two-body potentials, the interaction geometry is simply the square

of the distance between pairs of central atoms and neighbors. The interaction force is

found from the force table, which is subsequently broken into components and summed

on each particle.

When the three-body component of the force function is applied, the interaction

geometry becomes the distances between the central atom and two of its coordinating

neighbors of a certain type, and the angle between the two "bonds." Since the number

of interaction variables is large, but the actual number of interactions is rather small, a

force table for a three body potential would be difficult to compile, but fortunately, is

not required. The force is calculated for all three body interactions as they are needed,

and components are summed on each interacting particle.

As a particle moves, it may leave the primary cell. When this occurs, the particle

is translated back into the primary cell by subtracting the appropriate cell dimension and

adjusting for shear, as required. Thus particles that leave through one side of the cell

reappear through the opposite side of the cell, perhaps shifted if the cell has undergone

shear. This allows one to control the dimensions of the cell while permitting the atoms

to move as their environments dictate.

Observing the System

As a simulation progresses, one can collect information to find out what the

system is doing. Any aspect of a simulation can be recorded. A few that are commonly

recorded are instantaneous or average thermodynamic properties, particle positions and

velocities, pair correlation functions, and bond angle distribution functions.

System thermodynamics

Thermodynamic properties, such as system temperature, energy, and pressure are

calculated. A running average of these values may be kept, so that outputs of either

instantaneous thermodynamic properties or average properties may be made.

The temperature of the system can be found by assuming that the system behaves

as a classical statistical system, with average particle kinetic energies given by the

equipartition theorem, EK = k T. The temperature can then be adjusted by scaling the

velocities of all particles to give the proper average kinetic energy.

The system pressure or stress tensor is calculated using the Virial of Clausius,

which has the form

P +(NknT rI ENrVw(ro])

When converted for use in simulations, it takes the form

P O 3 NN