EXPERIMENTAL AND MOLECULAR DYNAMICS DETERMINATION OF

FRACTAL FRACTURE IN SINGLE CRYSTAL SILICON

BY

YUEH-LONG TSAI

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

1993

ACKNOWLEDGMENTS

I would like to express my sincerest gratitude to Dr. J. J. Mecholsky, my

supervisory committee chairman, for his invaluable guidance, support and encouragement

throughout the research portion of this project, and for his numerous suggestions and

corrections during the preparation of this dissertation. I am also thankful to the rest of my

committee members, Dr. R. T. DeHoff, Dr. F. Ebrahimi, Dr. P. A Mataga, and Dr. J. H.

Simmons, for their encouragement and advice throughout the research of this project.

I would like to express my special thanks to T. P. Swiler, from whom I learned a

lot about the MD technique and UNIX computer operation system. He also allows me to

use his computer codes to operate the MD simulation and MD movies.

I would also like to extend my thanks to Dr. T. J. Mackin, Mr. Z. Chen, Mr. L.

Hehn, Mr. J. Niaouris and A. Naman, with whom I had many precious discussions on

various aspects of this research. Special thanks must also be given to Ms. J. Y. Chan for

her help in obtaining some of the experimental data.

Finally, I am deeply indebted to my parents for their love and encouragement

throughout my entire life and to my wife, Yu, for her patience, understanding, and

support during the course of this work. I am especially grateful at this time because my

daughter, Gina, arrived in the world and she makes my school life more fruitful.

TABLE OF CONTENTS

Pages

ACKNOWLEDGMENTS......................................................................... ii

A B STR A C T ................................... ...... ...... ................ ..... .................... v

CHAPTER

1. IN TR O D U C TIO N .............................................................................. 1

2. FUNDAMENTAL.............................................................................. 7

Structure of Single Crystal Silicon ..................................................... 7

Failure A analysis ............ ....................... ........................ .. .............. 9

Fracture M echanics .............................................................. 9

Fracture Surface Analysis....................................................... 11

Indentation Fracture Mechanics......................................................... 15

Fractography-Indentation Analysis............................................ 15

Strength-Indentation Technique................................................ 18

Fractal G eom etry .......................................................................... 19

Fractal D im ension ................................................................ 19

Fracture Surface Analysis by Fractal Geometry............................ 27

Molecular Dynamics Simulation........................................................ 29

3. MOLECULAR DYNAMICS TECHNIQUES............................................. 33

O overview .................................................................................... 33

Potential Determination .................................................................. 34

Fundamentals About the MD Simulation.............................................. 36

Initial Conditions ................................................................. 36

Interactions Between Atoms.................................................... 39

MD Simulation Procedures .............................................................. 41

The Concept of Time Step...................................................... 42

The Method to Update Particle Positions.................................... 42

Periodic Boundary Condition .................................................. 43

Interactions Due to Two-Body Potential ..................................... 44

Interactions Due to Three-Body Potential.................................... 44

Updating of the Atom Positions ............................................... 45

Verlet's Algorithm ...................................................... 45

Gear's Algorithm........................................................ 46

Temperature Calculation........................................................ 48

M .D M ovie .................................. ....... ..... .... .......... ........... 48

Potential Energy of the System ................................................ 49

System Pressure Calculation ................................................... 49

Pair Correlation Function....................................................... 50

Bond Angle Distribution ........................................................ 51

Introduction of Strain............................................................ 52

4. EXPERIMENTAL AND SIMULATION PROCEDURE .............................. 53

Experimental Procedure.................................................................. 53

Sample Preparation............................................................... 53

Fracture Surface Analysis....................................................... 55

Toughness Measurement........................................................ 55

Fractal Dimension Determination............................................. 57

M.D. Simulation Procedure ............................................................. 61

Determination of Input Data.................................................... 61

Simulation Procedure ............................................................ 63

Fractal Analysis Using Simulation Results .................................. 67

5. RESULTS AND DISCUSSION.............................................................. 71

Experimental Results ..................................................................... 71

Fracture Surface Analysis....................................................... 71

Fracture Toughness Measurement............................................. 90

Fractal A analysis .................................................................. 100

Relationship Between Toughness and Fractal Dimension................. 104

Results From Molecular Dynamics Simulation ...................................... 111

Validity of the Applied Potential .............................................. 111

Fracture Using MD Simulation................................................ 120

Toughness in Different Orientations.......................................... 128

Strength Dependence on the Strain Rate ..................................... 130

Strength Dependence on the Crack Size...................................... 134

Fractal Analysis Using MD Results........................................... 142

Comparison Between Measured And Simulated Results ........................... 146

Comparison of Toughness ..................................................... 148

Comparison Between Experimental and Simulation Results for

Fractal Dimension............................................................... 148

6. C O N C LU SIO N S ............................................................................... 151

APPENDIX A. SUMMARY OF DATA FROM DIFFERENT INVESTIGATORS 154

APPENDIX B. STATISTICAL ANALYSIS................................................. 155

APPENDIX C. YOUNG'S MODULUS CALCULATION............................... 158

R E FERE N C E S .................................................................................... 160

BIOGRAPHICAL SKETCH .................................................................... 164

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

EXPERIMENTAL AND MOLECULAR DYNAMICS DETERMINATION OF

FRACTAL FRACTURE IN SINGLE CRYSTAL SILICON

By

Yueh-Long Tsai

December 1993

Chairman: Dr. J. J. Mecholsky

Major Department: Materials Science and Engineering

Single crystal silicon was selected as a model material in which to study the

correlation of the fracture surface features as characterized by their fractal dimension for

different orientations of fracture with the fracture toughness of the material as measured

using the strength-indentation, and fracture surface analysis techniques. Single crystal Si

was selected for several reasons: Si is a brittle, monoatomic material that will obviate the

complication of microstructure, i.e., grains, pores, etc. And Si has been well studied so

that many of its properties are well characterized. Flexure bars were indented with a

Vickers indent at various loads and fractured. After calculating the fracture strength and

toughness, the surfaces were analyzed and characterized using fracture surface analysis

and slit-island analysis. These analyses provided the size of the fracture initiating defect,

the geometry of the surrounding topography including the location of the regions of crack

branching and the fractal dimensions of selected areas on the fracture surface.

The { 100} plane is found experimentally to be the one with the highest toughness

of the three studied planes. The fracture plane on the {100} and {110} fracture plane

both initiate on the original plane and have the tendency to deviate to the { 111 } plane.

The fracture surfaces of Si have been analyzed and are found to be fractal both at

the atomic scale using the scanning tunneling microscope in a previous study and at the

micrometer scale using the optical microscope in this study. The irregularity of the

fracture surface is too complicated for Euclidean geometry to describe easily. Self-

similarity and scale invariance is suggested by the fracture surface appearance. The

fractal dimension is higher for the fracture plane with the higher toughness value. It is in

agreement with other studies which found higher toughness values associated with higher

fractal dimensions.

We also demonstrate the formation of fractal surfaces during fracture using a

molecular dynamics (MD) approach in a single crystal silicon structure and compare the

simulated results with experimental work. MD simulations using Stillinger-Weber

potential and Coulombic (Modified Born-Mayer) potential are performed in different

orientations to investigate the difference in toughness and fractal dimension with respect

to orientations. The close agreement between simulated and measured fractal structures

of the fracture surface suggests that this is a promising method for investigating atomic-

level processes during fracture.

CHAPTER 1

INTRODUCTION

The strength of brittle materials is determined by the stress required to

spontaneously propagate a crack. The strength is dependent on the local stress

distribution around a crack tip which is often treated using a stress-intensity factor

approach.1 The fracture toughness, which is a measure of the resistance to crack

growth, defines the work or mechanical energy expended in propagating the crack.

The fracture toughness can be represented by the critical stress intensity factor, the

fracture energy or the critical strain energy release rate. Techniques to determine the

fracture toughness of ceramic materials include cleavage,2 double cantilever beam

(DCB),3 crack indentation,4 strength-indentation,5 and fracture surface analysis (FSA)6

techniques. The mechanism of fracture in brittle materials is strongly influenced by

processes at the atomic level. The subject of fracture in brittle materials has long been

of interest; however, the underlying processes have remained difficult to investigate

directly. Some factors such as phonon-phonon interactions7 have been analyzed

recently, but the dynamics of bond breaking in brittle materials is only beginning to be

studied.8-11 Although it is intuitively obvious that there must be a relationship between

the bond breaking process and the macroscopic measures of fracture such as fracture

toughness and fracture topography, there has been little quantitative evidence of the

connection. Fractal geometry, which is a non-Euclidean geometry exhibiting self-

similarity and scale invariance, is a new tool that can be used to relate atomic scale

processes to macroscopic processes and features. The fracture surfaces of Si have been

analyzed and are found to be fractal both at the atomic scale12 using the scanning

tunneling microscope and at the micrometer scale13 using the optical microscope.

Fracture markings on glasses and polycrystalline brittle materials, known as

mirror, mist, and hackle, are precursors to crack branching14 and can be used to

describe the stress state6 and characteristics of crack propagation.15 These marking

have been observed for more than 50 years and were related quantitatively to the stress

condition in the 1950's.16 More recently, the repetition of these features was observed

and quantitatively related to stress intensity. Ravi-Chandar and Knauss14 also noted

that mist and hackle are self-similar; i.e., they appear to be physically similar and

produced in the same fashion. Their description, however, did not emphasize the self-

similar nature of the features.17 The distance to branching is directly related to the size

of the fracture initiating crack. Thus, one can imagine that the fracture surface features

can be related to the bond breaking processes at fractures.

Mandelbrot18 synthesized a branch of mathematics that provides a tool with

which to analyze self-similar and self-affmine processes. A self-similar process is one in

which a feature at one magnification is related to another at another magnification by a

scalar quantity. A self-affine process is one in which this magnification factor is a

vector quantity. Mandelbrot called this branch of mathematics "fractal geometry."

Mecholsky et al. 17 summarized that fracture in brittle materials can be described

as a self-similar process and thus, can be mathematically described by fractal geometry.

They derived an expression which relates the fracture of bonds at the atomic level (y) to

a characterization of the resulting morphology on the fracture surface of polycrystalline

materials, D*, which is the fractional part of the fractal dimension. They also derived

a relationship between the surface energy and the fracture energy of single crystals,

polycrystalline ceramics, and glass-ceramics. Their results suggest that modeling

brittle fracture as a fractal process may be useful in distinguishing toughening

mechanisms and in relating atomic bonding and fracture-surface morphology. The

present research was undertaken to study the relationship between fracture surface

morphology of single crystals and their resistance to fracture.

Unlike mathematically generated fractals, the self-similarity of fracture surface

geometry is bounded by measurement limitations. These limitations will be discussed

later. However, the self-similar nature of fractal surfaces offers a means of scaling

macro-scopically analyzable structures to microscopic processes on the atomic scale.

In addition, the fractal dimension has been shown to relate directly to fracture

toughness, thus forming a link between average macroscopic behavior and possible

atomic processes. 13'17

Single crystal silicon was selected as a model material in which to study the

correlation of the fracture surface features as characterized by their fractal dimension

for different orientations of fracture with the fracture toughness of the material as

measured using the strength-indentation, and fracture surface analysis techniques.

Single crystal Si was selected for several reasons: Si is a brittle, monoatomic material

that will obviate the complication of microstructure, i.e. grains, pores, etc. And Si has

been well studied so that many of its properties are well characterized. Flexure bars

were indented with a Vickers indent at various loads and fractured. After calculating

the fracture strength and toughness, the surfaces were analyzed and characterized using

fracture surface analysis and slit-island analysis. These analyses provided the size of

the fracture initiating defect, the geometry of the surrounding topography including the

location of the regions of crack branching and the fractal dimensions of selected areas

on the fracture surface.

In this study, we demonstrate the formation of fractal surfaces during fracture

using a molecular dynamics (MD) approach in a single crystal silicon structure and try

to compare the simulated results with the experimental works. The close agreement

between simulated and measured fractal structures of the fracture surface suggests that

this is a promising method for investigating atomic-level processes during fracture.

The use of molecular dynamics for studying the failure of materials has recently

become of great interest to investigators. Works by Paskin et al.,19 Dienes and

Paskin,20 and other collaborative works involving Paskin, have performed simulations

of crack growth in two-dimensional 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, and potential parameters. These works provide an atomistic basis to fracture

mechanics; however, they do not provide details on the fracture of particular materials.

271

Soules and Busby21 used molecular dynamics to study the rheological properties

and fracture behavior of a sodium silicate glass. Among the experiments in this work,

a uniaxial tension was applied to a system and the resulting behavior was studied. The

periodic boundary conditions used, i.e., on only two sides, made this experiment

equivalent to the fracture of a 50A diameter fiber, so bulk properties were not studied,

nor even approximated. This work was followed by works of Soules,8 Ochoa and

Simmons, and Ochoa et al.,10 which demonstrated the importance of dynamic effects

in the fracture of solids. Works by Ochoa, et al., 10 made use of full periodic boundary

conditions, i.e. on six sides, thus approximating fracture within a bulk material. These

simulations collectively studied the overall response of a system to an applied strain,

but did not study the atom-by-atom process in detail. A work by Simmons et al.11

examined the mechanisms involved in an individual bond fracture to describe a process

by which an entire fracture surface may be stabilized. Thus, these studies give possible

details on the fracture process for silicate glasses, and infer details of fracture for other

amorphous solids.

The difference in fracture between amorphous silica and single crystal silicon is

that silica has no preferred cleavage plane while single crystal silicon prefers fracturing

on the {111} and {110} plane.22 Thus, different loading orientations under different

loading histories in single crystal silicon would result in different fracture

topographies13'22 (although the fracture processes may be the same for different

orientations). The effect of elastic anisotropy on fracture provided the inspiration to

use single crystal silicon for the molecular dynamics simulation. From the simulated

loading history, the comparison of fracture toughness in different orientations can be

made. From the simulated fracture topography, the irregularity of the fracture surface

in different orientations can be analyzed using fractal analysis and the fractal dimension

can be compared with the toughness. Kieffer and Angell23 used the MD simulation

method to generate stable silica aggregates at various low densities similar to those of

experimental aerogels. They found that fractal dimensions and range of self-similarity

can be extracted from the radial distribution functions in those structure obtained from

MD simulation. They believed that a non-integer dimension is a characteristic feature

of the aerogel structure used in their study. However, because the way which they

used the change of the slope of the radial distribution functions with respect to the

density to obtain fractal dimension in their study was not reasonable, thus they got a

weird result which shows that a denser aerogel has a lower fractal dimension. It is

contradict to the fact that in the nature a full occupied volume has a dimension 3 while

a volume with porosity has dimension less than 3. Thus we must use a different

approach from theirs to obtained the fractal dimension in this study. Here the fractal

dimension on the generated fracture surfaces at the atomic scale from MD simulation

will be analyzed using slit-island technique and Richardson-plot.

This investigation will address several topics. First, there will be an

experimental determination of the fracture toughness in single crystal silicon as a

function of orientation of the crystal plane. Second, the fractal dimension of the

fracture surfaces of the single crystal silicon for the different orientations will be

determined experimentally. The fractal dimensions will be determined using

Richardson plots and slit-island analyses of contours of (horizontal) sections through

6

the fracture surfaces. Third, MD simulations using Stillinger-Weber24 and Coulombic

potentials will be performed in different orientations to investigate the difference in

toughness and fractal dimension with respect to orientations. Finally, there will be

discussions about the relationship between the fracture toughness and the fractal

dimensions of single crystal silicon for the different orientations and also comparisons

between experimental work and simulation results.

CHAPTER 2

FUNDAMENTALS

Structure of Single Crystal Silicon

Silicon is bonded by four covalent bonds which produce a tetrahedron. Each

silicon atom is bonded with four other silicon atoms due to the nature of the covalent

bonding. The covalent bond formed by two Si atoms sharing electrons is very

localized and directional. The space lattice of the diamond structure is face-centered

cubic with two atoms per lattice site, one at 0,0,0 and the other at 1/4, 1/4, 1/4 (Fig.

2.1). Each atom is tetrahedrally bonded to four nearest-neighbors due to the nature of

covalent bonding. As these tetrahedral groups are combined, a large cube is

constructed, which is called a diamond cubic (DC) unit cell. The large cube contains

eight smaller cubes that are the size of the tetrahedral cube but only four of the cubes

contain tetrahedrons. The lattice is a special face-centered cubic (FCC) structure,

which is shown in Fig. 2.1. The atoms on the corners of the tetrahedral cubes provide

atoms at each of the regular FCC lattice points. Four additional atoms are present

within the DC unit cell from the atoms in the center of the tetrahedral cubes.

Therefore, there are eight atoms per unit cell. The lattice parameter length is 5.43 A

and the unit cell diagonal length is 9.41 A as shown in Fig. 2.1. The Si bond distance

is 2.34 A.

Single crystal Si is generally classified as a brittle material whose atoms are

strongly covalently bonded. The primary cleavage plane is the {111} plane. The

principal factor for crack initiation is bond rupture.

(c)

Fig. 2.1 (a) Tetrahedron and (b) the diamond cubic (DC) unit cell. (c) The length

of AB is 5.43 A, AC is 9.41 A, and AD is 7.68 A.

J- I -- I --

Fracture of single crystal Si on the {100}, {110} and {111} plane has been

observed.13'22 Not only the toughness, but also the fracture surface topography were

found to be different on the different planes. These observations mean that elastic

anisotropy plays an important role in fracture; the bond breaking process is considered

the same since only one kind of Si-Si covalent bond exists.

Failure Analysis

The characterization of flaws and their relationship to strength plays an

important roll in the analysis of brittle materials. By comparing fracture mechanics

relationships with observed fracture surface markings, important information about the

fracture toughness can be determined. This section will present the basic principles of

fracture mechanics and fracture surface analysis of brittle materials.

Fracture Mechanics

One of the first analysis of fracture behavior of components that contain sharp

discontinuities was developed by Griffith.25'26 This analysis was based on the

assumption that incipient fracture in ideally brittle materials occurs when the magnitude

of the elastic energy supplied at the crack tip during an incremental increase in crack

length is equal to or greater than the magnitude of the elastic energy at the crack tip

during an incremental increase in crack length. By using a stress analysis developed by

Inglis, Griffith related fracture stress to the flaw size by means of an energy balance.

The relationship can be written as:

S=(2Ey1/2, (2.1)

whr -istC )

where arf is the failure stress,

E is the modulus of elasticity,

y is specific surface energy, and

c is the critical flaw size.

Although Griffith was the first to analyze the relationship between strength and

flaw size, it was Irwin who developed fracture mechanics into the present day form.

Irwin27-29 followed the work done by Griffith, Orowan,30 and Inglis31 to develop what

is known as linear elastic fracture mechanics. Irwin first analyzed the fracture of

flawed components using a stress analysis based on the Westergaard32 solution of an

elliptical crack in an infinite plate. For a surface cracked specimen under mode I

(tensile) loading, Irwin derived the following:

5f =, B (2.2)

where af is the failure stress,

Kic is the critical stress intensity factor for mode-I loading,

B is a geometrical factor which accounts for flaw shape, location,

and loading geometry, and

c is the critical flaw size.

Irwin also analyzed fracture in terms of the strain energy release rate, G. G is

defined as the elastic energy per unit crack length, U/c, and can be related to the

failure stress:27

= (Ecr' 1/2

CTf= -Ec (2.3)

where Gc is the elastic energy release rate at fracture.

The strain energy release rate can be related to surface energy, y, and it is as

follows for plane-stress condition:

G = 2y. (2.4)

In the above equation, only surface energy, y, is referred. Indeed during fracture

process, the released strain energy not only produce two free surface, but also produce

different terms like sound, light, electron emission. Thus, the y used here indeed is not

just as simple as surface energy and some other energy is also included.

The stress intensity factor can be related to the strain energy release rate and the

fracture surface energy. By combining those previous equations, it can be shown that

Kic= = EG= -2yE. (2.5)

Fracture Surface Analysis

Useful information on the mechanical and fracture behavior of a failed ceramic

can be gained from microscopic examination of the fracture surface.33-35 The fracture

surface of brittle materials exhibits distinct fracture markings surrounding the critical

flaw. These features shown in Figure 2.2 can be related by known fracture mechanics

relationships to provide additional information on the fracture process and can also be

used to locate the origin of failure.35

The fracture surface markings shown in Figure 2.2 can be used to describe the

stress state of a brittle material.14'33'35'36 Four different regions referred to as the

fracture mirror, mist, hackle and crack branching can be seen in Figure 2.2. These

regions are associated with particular stress intensity levels and crack velocities which

are responsible for each distinct region.34 As a crack propagates from the critical flaw,

a smooth region which is basically perpendicular to the tensile axis is formed. This

smooth region is known as the fracture mirror. The fracture mirror is typically

bounded by a region of small radial ridges known as mist and the mist is bounded by

an even rougher radial ridge region called hackle. Finally, the propagating crack

reaches a characteristic energy level and crack branching occurs. Crack branching is a

region where two or more cracks form from the primary crack front.

Fig. 2.2 Schematic of features found on the fracture surface of a brittle material

subjected to a constant load. Solid semi-elliptical line at center represents the initial

flaw size depth ai and width 2bi. Dashed line represents the outline of critical flaw

depth acr and width 2bcr. Mirror/mist, rl, mist/hackle, r2, and crack branching, r3,

radii are shown along the tensile axis.

r,2

--'cb7------ *

Once the failure origin is located after fracture, its size can be measured and

used to determined the toughness of the material. Figure 2.3 shows a semi-elliptical

crack in a plate. It can be shown through a Griffith/Irwin approach that the stress

intensity on a semi-elliptical surface crack of depth, a, and half-width, b, can be given

as:35,37

K, = 1-.--. _27 Y-aJa (2.6)

where a is the semi-minor axis for an elliptical crack,

ca is the applied stress, and

D is an elliptical integral of the second kind. D varies between 0 = 1

for a slit crack (a/b =0) to 0D = 1.57 for a semi-circular crack (a/b= 1).

The criterion for failure in a brittle material is that Ki _> Kic at which point the

initial flaw will propagate spontaneously, where KI, the stress intensity, is a measure of

the magnification of the external loading at the crack tip, and KIc is the critical stress

intensity factor or the fracture toughness of the material. The strength or failure stress,

cf, of a brittle material can be related to the flaw size:

a DK c (2.7)

For a semi-circular surface crack, which is stress free and which is small relative to the

thickness of the cross-section, ) = 1.57 and Kjc can be found as:

Kic=1.24 cyfVc, (2.8)

where c is the square root of a-b and is the size of a semi-elliptical crack which has the

size of an equivalent semi-circular crack.30

Thus, the fracture toughness of a material can be calculated if the failure stress

is determined and the crack dimensions a and b are measured.

14

a

CY

2b

Fig. 2.3 An semi-infinite plate under uniform tension and containing a semi-

elliptical crack.

Indentation Fracture Mechanics

Indentation fracture mechanics has been well established as an important method

for the study of the mechanical behavior of brittle materials.38-43 Indentation provides

a means of introducing flaws into a material with a controlled size, shape, and location

on the sample. This is in contrast to materials with naturally occurring flaws where the

flaws vary in size, shape, location and concentration. Indentation provides a means of

characterizing the fracture process by introducing a controlled flaw into a brittle

material. The following discussion will concentrate on fracture surface analysis and

strength indentation analysis.

Fractography-Indentation Analysis

Figure 2.4(a) shows a schematic of the indentation deformation/fracture pattern

for the Vickers diamond geometry: P is the peak load and a and b are characteristic

dimensions of the inelastic impression of the radial/median crack, respectively. Figure

2.4(b) shows the schematic of a fully developed Vickers indent. The flaw is taken to

possess a penny-like geometry. The following discusses the cases of zero and non zero

residual stress terms separately:40

(1) Zero residual contact stresses: If the indentation flaw were to be free of

residual stresses, the stress intensity factor for uniform tensile loading would have the

standard form

Ka = YfC12, (2.9)

where Ka is the applied stress intensity factor without residual stress,

Y is a unitless geometrical factor. For a semi-circular crack, Y = 1.24,

C

J

-1

I-

Fig. 2.4 Configuration of median/radial cracks for Vickers indentation showing:

(a normal indent load, P, generating median opening forces Pe (elastic component) and

Pr (residual component), (b) load removal eliminates the elastic component, and (c)

fully developed radial crack pattern.

cf is the failure strength, and

c is the critical flaw size.

(2) Nonzero residual contact stresses: With residual contact stresses present at

the flaw origin it becomes necessary to include an additional tensile term in the stress

intensity factor:44

KI = Ka + Kr = Ycac2 +XrP / c32, (2.10)

where Kr is the stress intensity factor with residual stress,

P is the peak load, and

Xr is a unitless geometry constant.

For large c values, the applied-stress term controls the fracture, as before, for

small c values, it is the residual-stress term that dominates. "Under equilibrium

fracture conditions the flaw will accordingly undergo a precursor stage of stable growth

as the tensile loading is applied; failure then occurs when the crack reaches a critical

size, at which point the applied stress is intense enough to cause spontaneous

propagation."40 This critical configuration is obtained by inserting K, = Kc (Kc is

used here instead of Kic because, strictly speaking, the toughness determined using

indentation is different than Kic; however, in practical terms, they are equivalent) into

the above equation and evaluating the instability condition dca/dc = 0:

= 4Xrp )2/3 =(3Kr_ )2 (.1

c Kc 4Ycrf(2.11)

from which Kc can be found as:45,46

4 1/2

Kc Yyfc2. (2.12)

3

For a semi-circular flaw, which is small relative to the thickness of the material into

which it is placed, Y = 1.24, and the above equation can be rewritten as:

Kcl = 1.65acfc12. (2.13)

Thus, the resistance to fracture can be determined from a measurement of the surface

crack size and the strength of bars that have been indented.

Strength-Indentation Technique5

Consider the Vickers diamond-induced radial cracks, which are shown in Fig

2.4. If the applied loading is uniaxial, the indentation is aligned with one set of

pyramidal edges parallel to the tensile axis.

It has been shown that the toughness may be derived using the following

expression:

-1/8 "/

V H)

K =nv T1{J (cafP1/'3) (2.14)

f')76 y23R1/4

where 11R = 256Y2/3 is a geometry constant,

E is the Young's modulus, and

H is the hardness.

Replacing rv by an average quantity, 0.59, would add no more than 10% to the error

in the Kc evaluation for a material whose elastic/plastic parameters are totally

unknown. So the above equation can be rewritten as

(gE1/8 )3/4

Kc = 0.59 E ((fP1/3) (2.15)

Thus, the resistance to fracture can also be measured from a flexural test in which the

specimens have been indented at different load levels without measuring the crack size.

Fractal Geometry

Fractal geometry is a non-Euclidean geometry that was rediscovered,

popularized and applied by B. B. Mandelbrot.18 The word fractal is derived from the

Latin "fractus" which means fragmented or broken. Nature exhibits diverse structures

having inherent irregularities. Some of these structures can be described using

Euclidean geometry but other structures are better described by using fractal geometry.

In his book, The Fractal Geometry of Nature,18 Mandelbrot describes the extensive

applications of fractal geometry. The concepts have been used to described the

geometry of clouds, soot aggregates, dielectric breakdown, coastlines, and other

natural phenomena, including fracture.

Fractal geometry exhibits self-similarity or self-affinity, shows scale invariance,

and is characterized by the fractal dimension. Self-similarity means that features in

different regions appear to be similar to one another. They can be related by a scalar

multiple. Features in a self-affine object are related by a vector quantity. Scale

invariance occurs if a feature on different scales appears to be the same. These

properties indicate that a fractal can be created by using a shape to be repeated and

reduced in size following a prescribed sequence.

Fractal Dimension

In Euclidean geometry, objects occupy integer dimensions, i.e., 1, 2, 3, and so

on. Fractal geometry admits to the description of objects that occupy fractional

dimensions. Fractal geometry is characterized by its fractal dimension which is related

to measurement theory.

Consider how we measure the length of a line. The measurement of a line

requires a measuring tool or ruler. A line is given in Fig. 2.5(a). The length of this

line is measured by covering it completely with as few of a given radius, R, discs as

possible, as in Fig. 2.5(b). The length is given by

length = (number of covering discs) x (2R). (2.16)

If we have a new supply of discs with a different radius, say R/2, then the line length is

determined in the same manner, and a new line length is computed, as shown in Fig.

2.5(c). Note the difference in the measured line length in Fig. 2.5(c). Although the

measuring disc was decreased to half of the original measure, the number required to

cover the line more than doubled. The "tortuous" nature of this curve gives a length

that is scale dependent. As the measuring scale decreases, the measured length

increases.17

The measured length of a line can be described by the following equation,

called Richardson's equation:

L = kSID, (2.17)

where L is the measured length,

S is the measured scale,

D is the dimension of the curve,

k is a proportionality constant.

If D is equal to one, the line is Euclidean and its measured length is not a

function of scale. If 1 < D < 2, the curve is said to be fractal with a dimension given by

D. Notice that for a given line segment, the length is not constant, but the dimension

is. Thus, the fractal dimension is the dimension in which a measure can be made, not

the measure itself.

The fractal dimension of a line can be computed from the above equation. The

line's length is computed over a range of scales. A log-log plot of length vs. scale

gives a straight line with slope equal to 1-D, as shown in Fig. 2.6.

(a)

O diameter =2R

Fig. 2.5 Measurement of a line segment ( or signal). (a) A schematic for a line

which has texture. Measuring the length of a line requires a ruler. The (circle)

covering has a diameter, 2R. (b) The line is completely covered to determine its

length. (c) Changing the size of the measuring disc may produce a different measure of

the length.

Continued.

Q diameter = 2R

length = # of disc x diameter 2R

O diameter = 2r

length -= # of disc x diameter 2r

Fig. 2.5

Slope = 1-D

log L

LogE

Fig. 2.6 A log-log plot of length vs. scale gives a straight line with slope equal to

1-D.

Mathematical fractals will obey the above equation for all scales. However,

physical fractals may have a finite cutoff at both large and small scales. In the physical

world, the fluctuation of a line will find its limit in the smallest measurable feature.

The ruler length determines the smallest feature we can measure. If variations in

structure are below the sensitivity of the measurement scale, the line will appear to be

Euclidean. In that case, the longest line we could measure would depend on the finest

possible scale of measurement. We could not conclude that the line has no

fluctuations, only that we have no ability to measure the fluctuations.

On the contrary, a line may be composed of indivisible components. After our

measuring stick has become smaller than the smallest of these components, the line

would cease to exhibit scale-dependent length. This length, then, would represent a

true upper bound. Such a line could be called "fractal" over the range, it exhibited

scale-dependent length. It is not, however, a fractal in the strict mathematical sense.

The combination of the two limits would result in the Richardson plot as shown

in Fig. 2.7. The object would be described as a fractal over a bounded range of scales.

The cut-off scales are of special significance to the particular geometry of a specific

physical phenomenon.

A box is shown in Fig. 2.8(a)18'47 whose sides are of unit length. Next to the

box is a shape called a "generator," composed of line segments scaled to 1/4 of the

length of a side of the box. Each side of the square is replaced with this "scaled"

generator. The result is shown in Fig. 2.8(b). Again, the generator is scaled down to

the length of each straight line segment of this new object. The straight line segments

of Fig. 2.8(b) are replaced with the scaled generator. A portion of the resulting shape

is Fig 2.8(c). This process is continued, generating an object which is said to be scale

invariant and self-similar. A scaling fractal looks geometrically the same everywhere

and on all scales.

Fractal on bounded

range

log L

log E

Physical fractals have scale dependence over a bounded range of scales.

Fig. 2.7

The stepwise construction of a scaling fractal.

each side of the square

is replaced by the scaled

generator.

again, the generator is

scaled to replace each

of the segments above.

h-j

Fig. 2.8

The fractal dimension of this object, however, can be computed in another,

quite different fashion. Scaling fractals obey the following:

NrD = 1, (2.18)

where N is the number of elements in the generator,

r is the scale factor of an element,

D is the similarity or fractal dimension.

For the plot of Fig. 2.8 called the quadratic Koch curve, N=8, r= 1/4, Thus,

the fractal dimension D = log8/log4 = 3/2.

Fracture Surface Analysis by Fractal Geometry

Mecholsky et al.17'48 experimentally determined a relationship between fracture

toughness (Kc) and fractal dimension (D). They found that

Kc =A (D -1)1/2 =A D*1/2, (2.19)

where D* is the fractional part of the fractal dimension,

Kc = Kic Ko, where Ko is the toughness of the material for a smooth

(Euclidean) fracture surface, and

A is a constant.

Thus, as fractal dimension increases, fracture toughness increases. In some

papers (D-l) is referred to as D*. Thus the above equation can be rewritten as

Kc = AD*1/2, (2.20)

where A is a family parameter. The family parameter, A, identifies a line in the Kc

vs. D* plane. Within a family, an increase in fractal dimension corresponds to an

increase in fracture toughness.

Since D* is dimensionless, a dimensional analysis of this equation requires that

A have the dimensions of toughness. Mecholsky et. al.17 proposed that A is a product

of a characteristic length and Young's modulus, so that A = E(ao)11/2. Then the above

equation can be written as:

Kc = E ( ao D* )1/2, (2.21)

where Kc is the fracture toughness,

E is the Young's modulus,

D* is the fractal dimension, and

ao is the characteristic length.

Recalling that Kc = EG = -2yE for a plane stress condition, ao can be

obtained by the following equation,17

ao =- 2y (2.22)

ED*

In order to prevent getting an infinite values of ao calculated from the above

equation if D* is approaching zero. A modified equation from a different approach to

find ao is expressed as:

Y = o70 + aoEhkD*/2, (2.23)

where y is the fracture energy,

7o is the surface energy for an Euclidean fracture surface (D* = 0),

ao is a characteristic (atomic) length, and

Ehk is the Young's Modulus of the hkl fracture plane.

Using experimentally obtained values of y, E and D*, yo and ao can be

calculated. The exact interpretation is not clear at this time; however, it is reasonable

to assume that it is related to the characteristic length of the generator. In terms of

fracture, this would correspond to a characteristic bond length which is involved with

initial fracture at the crack tip.

Mecholsky et al.6 reported that the outer mirror constant, i.e., oa(r2)112 where r2

is as noted in Fig. 2.2, is related to Kc and independently to E. Kirchner49 shows a

similar relationship with the crack-branching constant and used this information to

support his idea of a constant strain intensity at branching. Tsai and Mecholsky50 used

a fracture energy concept to prove that the mirror-hackle boundary forms at a constant

energy value. Mecholsky51 points out that the rj/c = constant. Freiman et al.52

showed that the relationship could be represented as:

Kbj= (bj)1/2E, (2.24)

where bj is a characteristic dimension on the atomic scale, and

j = 1, 2, 3 is related to mirror-mist, mist-hackle and crack branching

boundary.

They further showed that if b1 = ao, where b1 corresponds to the mirror-mist

boundary, then the flaw-to-mirror-size ratio is equal to the fractal dimensional

increment:

c/rI = D*. (2.25)

This equation can provide a link between fractography and fractal geometry.

They found equation (2.25) in good agreement with experimental results. Thus, it

appears that the relationship between the crack size and the formation of the branching

boundaries is related to the dimensionality of the structure. Presumably, the

dimensionality of the structure is related to the bonding, i.e., the strength and length of

bonds on the atomic scale.

Molecular Dynamics Simulation

Molecular dynamics simulation is a computer technique to model material

structure at the atomic level in order to model or study the material properties with an

atomic view point.53'54 During the molecular dynamics simulation, a small box with a

limited number of atoms is generated. The small box is considered a simplified system

which is ideally a realistic subset of the real system. This MD technique uses

Newton's equations of motion at constant acceleration over very short time intervals.

Interaction forces generated from the potential energy generated between the atoms are

assumed. The acceleration of particles (or atoms) results from the interaction of those

interatomic forces. The interaction forces produce the fundamental momentum for the

particles inside the small and idealized system.

Molecular dynamics simulations can be applied to study the properties of

materials. During MD simulations, the position and velocity of each atom can be

calculated to show the dynamic properties at that instant of time. The power of the

MD simulation is that the rule of interaction can be changed at will to model any

complex sets of particles, and environments can be simulated at conditions which are

not accessible experimentally. The difficulty with the MD simulation is that the rule of

interaction, i.e., the selected potential, must be carefully chosen in order to be realistic

or meaningful.

Three major uses of MD simulation are the development of theoretical models

for materials, the testing of theoretical models, and the predictions of material behavior

under experimentally inaccessible or costly conditions.

Material models can be developed from molecular dynamic simulations to

observe what occurs during the simulation and be able to set the rules of a simulation.

If we want to model the fracture of silicon, we do this by observing individual atomic

motions as samples are fractured using identical interaction potentials. We can also

fracture at several strain rates to investigate the effects of allowing the material to relax

by thermal vibrations as compared to not allowing the material to relax. After

performing such simulations, we can develop a general idea of how single crystal

silicon fractures and those ideas will be valuable to those studying the strength of Si.

Thus, molecular dynamics simulations can be used to devise models.

Molecular dynamics can be also used to test theoretical models using the ability

to observe all modeled atomic motions and the ability to simulate a wide range of

environments. If the system behaves as predicted by a model under a range of

conditions, then the model can get greater credibility.

Extreme conditions also can be simulated using molecular dynamics. As one

might simulate any environment, one needs only to have a realistic simulation to

predict materials behavior or study possible modes of failure in extreme environments.

Thus, as molecular dynamics becomes better at simulating reality, one may expect this

application of molecular dynamics to find regular use.

When creating a molecular dynamics system, one would try to model the laws

of nature as close to reality as possible in order that the processes observed in the MD

simulation can be as realistic as possible, and that the resulting character of the

simulated system corresponds to the character of real materials. Because the laws of an

MD system are executed by high speed computer of limited computational power, there

must be limitations of the resulting MD system. These limitations are seen as

limitations in the spatial and temporal resolution in nature, as well as limitations on the

number of interactions possible within such a system.

The spatial dimensions of a molecular dynamics cell must be within a certain

order of magnitude of the computer precision of positions obtainable.53 The spatial

resolution of the computer calculation is generally 1 in 107 or 1 in 1016, depending on

whether single or double precision is used in the computer calculations. Over long

calculation periods, the use of single precision has resulted in a small error in the

system, while yielding a lot of time-saving is obtained. Thus, here the choice of spatial

resolution is not a problem in performing the MD simulations.

The length of time represented by the simulation must be within a certain order

of magnitude of the time step, the basic unit of time in a simulation. Generally the

duration of an experiment is between 10,000 to 100,000 time steps depending on

computational resources and the purpose of the experiment. The temporal resolution of

an MD simulation may be a serious limitation. The reason for this is that processes

relevant to the structural evolution of materials have characteristic times spanning about

19 orders of magnitude, far more than the four to five orders of magnitude spanned by

MD simulations. Generally a simulation which lasts for hundreds or thousands of

atomic thermal vibrations is needed in order to understand the properties of the system.

The mass of the electron is about one ten thousandth of those of atoms. The

motion of electrons cannot be modeled with the motions of atoms at the same time

since the characteristic time is so small compared with that of atoms. Thus electronic

motions cannot be used to provide true interatomic potentials of materials in such a

simulation and the motions studied here are limited to those of atoms only. The chosen

potential is necessary to accurately represent the thermal vibration motions of atoms in

order to model the properties of the simulated system. Thus the smallest time step used

must be a small fraction of the vibration period. A typical time step chosen in this

study is 0.5 x 10-15 sec. One generally can only run simulations for about 105 time

steps which results in a limitation in the duration of modeled experiments to hundreds

or thousands of atomic vibrations.

Another limitation is that the number of interactions between atoms cannot be

too big. The more atoms put in the system, the more computer CPU time will be

consumed. For example, if a system with n particles is created and only one-to-one

interaction is considered, then the total interactions between those atoms is n2. The

time to calculate a system with 10On atoms is 102 longer than that needed to calculate

one with n atoms. Thus the simulated system must be kept small and the chosen

interacted potential must be simple and easily applied.

Thus the differences between an MD system and the real universe are that the

length of studied time and the size of the studied system. The MD system has to be

small because of the limits on the spatial resolution of time and number of atoms used

in the system.

CHAPTER 3

MOLECULAR DYNAMICS TECHNIQUES

Overview

A technique for studying materials whose popularity and usefulness has arisen

within the past few decades is molecular dynamics (MD) computer simulation. This

technique can solve equations of motion for a system of particles which may be ions

(ionic bonding) or atoms covalentt bonding or metallic bonding). Each particle moves

according to the forces on it caused by all of the other particles in the system. The

inter-particle forces are derived from an assumed potential function which describes the

interactions between the atoms. Several types of systems have been simulated by MD

with various interatomic potentials. Those systems include hard sphere systems, soft

sphere systems, ionic systems and Lennard-Jones systems. These systems are

generally simulated using simple pair interaction potentials to obtain a close

approximation.

The covalent bond formed by two Si atoms sharing electrons is very localized

and directional. In ionic systems, Coulomb potentials which include short range

repulsive forces as well as attractive forces, are required. In covalent systems, bond-

directionality should be included. Two body interatomic potentials are not enough to

produce the diamond cubic structure. Since Si exhibits that structure as a solid, three

body potentials must be introduced in order to conduct a reasonable simulation. The

Stillinger-Weber potential24 has been used for Si solid simulations using MD and the

simulated pair correlation function agrees very well with the experimental data. It was

also generally found able to accurately simulate the elastic and thermodynamic

properties of Si.55'56 This potential has been adopted here.

Potential Determination

Single crystal silicon consists of atoms held in place (a diamond structure which

is relatively anisotropic from the atomic viewpoint) by strong and directional bonds. It

seems reasonable at first sight that the corresponding potential (D could be

approximated by a combination of pair and triplet potentials, (2 and cD3. And these

two potentials can be expressed as a function of energy, F, and length units, 5. Thus,

02(rij) = s f2(rij/5), (3.1)

and

't3(ri, rj, rk) = s f3(ri/5, rj/5, rk/6). (3.2)

We use a reduced pair potential for Si selected from the following five-

parameter family, a simplified Lennard-Jones potential, from the work of Stillinger and

Weber:24

f2(r) = A (B r-P r-q)exp[(r a)-], if r < a; (3.3)

= 0, if r>a.

Thus,

(D2(rij) = e f2(rij/5) (3.4)

= A (B ()-P (r)-q) exp[( r )-I], ifr

6 6 6 6

= 0, ifr>a.

This form automatically cuts off at r = a without discontinuities in any r derivatives,

which is a distinct advantage in any MD simulation application. Stillinger and Weber

use the same cutoff technique for the three-body interactions, f3(ri,rj,rk) = hjik + hijk

+ hkij. The function h is given :

h(rij, rik, 0jik) = ) exp[ri /(rij a) + ril /(rik a)] x (cosOjik + 1/3)2, if r < a;

(3.5)

= 0, ifr>a.

where 0jik is the angle between rj and rk. The function h exists when both rij and rik

are less than the previously introduced cutoff a.

The parameters which best satisfy the radial distribution function were

determined by Stillinger and Weber to be: A = 7.049556277, B = 0.6022245584, p =

4, q = 0, a = 1.8, 1 = 21.0, h = 1.20, d = 0.20951 nm, and s = 50 kcal/mole =

3.4723x10-12 erg/atom pair.

Although the Stillinger-Weber potential has been successfully used to model the

structure of silicon, it does not provide the expected repulsive forces following covalent

bond breakage which result from ion formation during the fracture process. The Si-Si

covalent bond is formed by an attractive force even if the distance between the two

atoms is larger than the equilibrium position as long as it is within the cutoff distance.

However, during fracture this condition is not realistic and makes the Si-Si bond

difficult to break. An unrealistically long fiber-like structure is observed to result if the

unmodified Stillinger-Weber potential is used. An improvement developed here

assumes that following bond breakage each Si atom will gain one positive charge for

each lost neighbor. Si will have two positive charges if it loses two of its four

neighbors and four positive charges if all four neighbors are displaced beyond the

cutoff distance.

A modified Born-Mayer potential1'57 is used for modeling the resulting ionic

forces. Consequently this approach assumes ionic repulsive forces when two Si atoms

both with positive charges fall inside the Born-Mayer cutoff distance. This potential

has the following form:

-r qiqj r

rij = Aij.exp() + ( ). erfc(-), (3.6)

p r

where erfc() is the complement error function,

E is a constant,

L is a linear dimension of the molecular dynamics cell,

q is the charge,

r is the distance between the i and j atoms, and

p is a hardness parameter.

After combining the Stillinger-Weber potential and the modified Born-Mayer

potential, the fiber-like structure disappears and potentially more realistic fracture

surfaces can be obtained.

Fundamentals About the MD Simulation

How the MD simulations are performed will be explained now. First, we

discuss how the system was set up, which includes the initial particle positions,

velocities and the potential described in the previous section. Then we discuss how the

system is operated, which includes the application of periodic boundary conditions, the

force summation on the particles, control of the velocities of particles, and control of

thermodynamic values. Finally the discussion of how the data from the simulation

system is collected will be given. A flow chart in Fig. 3.1 is given to illustrate how

the MD simulation is operated.

Initial Conditions

Initial particle positions are located at specific positions since a single crystal

structure is desired.

Initial Velocity:

temperature

kinetic

k energy

Ek

Velocity (V)

random

distribution

function

Interactions between atoms:

apply strain

+(V+ AV) x At

t

new position

of atoms

force

between

atoms

I

acceleration D- AV

ax At

A flow chart for the MD simulation.

position

of atoms

Fig. 3.1

38

Initial velocities are given using the Poisson velocity distribution. The velocity

distribution is generally selected to correspond to a temperature at which the crystal is

stable. Here it corresponds to room temperature or to the temperature of the

simulation.

Random velocities are assigned to particles using the Poisson velocity

distribution. The formula giving the Maxwellian velocity distribution in three

dimensions are58

Vi 2kBT

Vi = 2kBT x ln(Randl) x cos(27' x Rand3) (3.7)

mi

I- M\i I

Vi= 2kBT x ln(Randl) x sin(27Tc x Rand3) (3.8)

Smi

V I- 2kBT

Vi= B x ln(Rand2) x sin(27 x Rand4) (3.9)

mi

where Vxi, Vyi, and Vzj are the components of the velocities of particle i,

Randl, Rand2, Rand3, and Rand4 are random numbers evenly

distributed between 0 and 1,

mi is the mass of particles i,

T is the temperature, and

kg is the Boltzmann's constant.

Particle velocities are applied by putting the velocities obtained above into

positions before and after the current time step, which is the basic unit of time during a

simulation. Thus the random velocities are effectively converted into average velocities

during the previous time step, and particles positions and velocities are given by the

difference of two sets of positions for each particles. The reason for this will become

clear in the discussion of Verlet's algorithm.

Interactions Between Atoms

Interaction forces are tabulated for all types of interactions possible over a range

of interparticle spacing. The interatomic potential assumed in this study is the

Stillinger-Weber potential and the modified Born-Mayer potential as illustrated above.

The derivative of these potentials is the force as a function of r. The two-body

potential is

2(rij) = 6 f2(rij/8)

=EA(B (6)-p ()q)exp[( 6 )-], ifr

= 0, ifr>a.

Thus, the two-body force is

f2(rij) = 2 (3.10)

ar

exp( )-1Y

8 r-a

(Bp (r)-p-1- (r)-q-1)+B(r)-P_( y 82)2

5 8 8 8 o (r -a)2J

,ifr
= 0, ifr>a.

It is computer-time-effective to tabulate values of this function as a function of

r2 and divide by r so that one may use the square of the interaction distance in finding

interparticle forces in the table. One may find the Cartesian components of the force

between two particles by simply multiplying the appropriate tabulated value by the

difference in Cartesian coordinates of the two particles.

40

The energy of potential energies for interactions between particles are also

tabulated in terms of r2. This is used for the calculation of potential energies for the

system.

The three-body force is considered in the same way.

D3(ri, rj, rk) = s f3(ri/5 rj/5 rk/5).

f3(ri,rj,rk) = hjik + hijk + hkij.

The function h is given :

hjik = h(rij, rik, 0jik) (3.11)

= X exp[h /(rij a) + Tj /(rik a)] x (cosOjik + 1/3)2, if r
= 0, ifr>a.

The following derivation related to the three-body potential will only use hjik for

convenience. Thus the three body force related to hjik is as follows:

a -, 1 a i

VO = gradD = er + 0 -- (3.12)

aT r aO

F=-V =-(e +e1 (3.13)

(A) At first, we differentiate r1 only.

F = Fr + F0 =-V(D (3.14)

At dr direction,

Fr- = --)2 exp[7y(rl -a)- +y7(r2 a)-l](coso+ )2 (3.15)

ar (r, -a) 3

At O direction,

F_=- 2X exp[7 (r1 a)- + y (r2-a)-1]sinO(cosO+-) (3.16)

r1 ao r1 3

Now transform the polarized coordinates into Cartesian coordinates.

Fr +Fe=Fx +Fy +Fz (3.17)

41

Consider Fr first. At first we have to find the transformation of u1 into Cartesian

coordinate.

U = UIx + uyj + ulzk (3.18)

FrUil = Fr(uixiJ+uiyj+uizk) (3.19)

Now consider F0. According to Fig. 3-1,

60 = fui x i3 = i1 x(i1 x i2) (3.20)

and because A x ( B x C) = (A.C)xB-(A.B)xC, thus

e0 =i x(il xu 2) (3.21)

=(i1.U2)x ii1 -(1il-i1)U2

=cosOUi1 iU2

Fo0e=Fo(cosOi1 -Ui2) (3.22)

= F(cos0uix U2x)i + FO(cos0uiy U2y)j + FO(cos0ulz -U2z)

Thus,

Fr +F0=Fx +Fy +Fz (3.23)

[FrUix + F0 (cosOulx- U2x)]1

+[FrUiy + F0(cos0Uly U2y)]j

+[Fruiz + F0 (cos Ouiz- U2z)]k

(B) We can consider the r2 term in the same way as (A).

MD Simulation Procedures

The explanation of how the simulations are performed will be given here.53

The concept of time step, the interaction between particles, the method used to update

particle positions, the application of periodical boundary conditions, and the control of

system temperature will be discussed as follows.

The Concept of Time Step

The time step is the basic unit of time during a molecular dynamics simulation.

As there are many atoms interacting in the simulated system, the motions of particles

cannot be solved analytically. Thus one can calculate the motions of atoms using

Newton's laws of motion where a time step dt is used in which the acceleration is kept

constant. A careful choice of time step is very important. If the time step is chosen to

be too small, a waste of computer running time results. If the time step is chosen to be

to big, nonlinear effects have to be considered and the atoms will move too far away to

reach unrealistic positions before the next time step is calculated. Following the work

of Soules, Ochoa, and Swiler, the time step used in our simulation is 0.5 x 10-15

second, and is about 1/40 of an atomic vibration period. A test to decide if a time step

is of sufficient length is to perform a simulation over about a thousand steps and check

to see if the system energy remains fairly constant. If not, the time step is too long. A

time step of 0.5 x 10-15 seconds was found to be of acceptable length and was applied

here.

The Method to Update Particle Positions

It is very important to handle particle interactions in a time-saving way.

Particle interactions are allowed only when particles are within a certain distance for

others. Beyond the distance, the effect of interaction is either none or negligible.

Cutoff distance is the distance in which interactions are allowed. Thus a neighbors-list

for all particles which lists all the neighbors which may be considered for the next 5 or

10 time steps will be made and updated at every 5 or 10 time steps depending on the

pulling rate. The fixed distance for the neighbors-list is chosen to be the sum of the

cutoff distance and the distance a particle may move in the next 5 or 10 time steps.

Only the neighbors listed on the neighbors-list will be considered when calculating

interactions.

The squares of the distances between particles are calculated using the

Pythagorean theorem for time-saving reasons. Square roots of these values are not

taken when calculating interaction forces because square root operations are time-

consuming and the force table is tabulated in terms of distance squared. Square roots

of the distances squared are only required when the neighbors-list are updated and

when pair correlation functions and bond angle distribution functions are accumulated.

Periodic Boundary Condition

Periodic boundary conditions are applied when finding inter-particle distances.

If periodic boundary conditions are applied, particles which leave the primary cell after

a time step will be translated back into it by subtracting the appropriate cell dimension.

Thus particles which leave through one side of the system cell will appear through the

opposite side of the cell.

Thus if the difference in the coordinate is greater than one-half the

corresponding dimension of the system cell, the corresponding image of the particle

located in the 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 which is used in the Pythagorean

relation to find the square of distance between the particles. The algorithm by which

this is performed is as follows:

dijk = Xij Xik (3.24)

dijk = dijk si trunc(2 ) (3.25)

r2 ,2 d'2 +,2

2 =dljk + djk (3.26)

rjk Ij 2jk + 3jk

where i = 1, 2, 3 denote the Cartesian dimension of space,

Xij is the i Cartesian coordinate of atom j,

Xk is the i Cartesian coordinate of atom k,

dijk is the i coordinate difference between atoms j and k,

dijk is the i coordinate difference between atoms j and k after correction

for periodic boundary conditions,

si is the i dimension of the system cell,

trunc() is a function which returns the integer closest to the argument

between the argument and zero, and

rij is the distance between atoms i and j.

Interactions Due to Two-Body Potential

Two-body interactions between particles are handled in a four step process.53

The square of the interparticle distance is found for particles located in the neighbors

lists. The interparticle force is then found by looking up the appropriate value in the

force table. The components of the force are found by multiplying the value tabulated

in the force table by the difference in the coordinates for each dimension (dijk). The

components of force are then summed for each particle for all interactions.

Interactions Due to Three-Body Potential

Three-body interactions between particles are handled in a four step process.

The square of the interparticle distance is found for particles located in the neighbors

lists. Two atoms, j and k, within the cutoff potential range are found from the

neighbor list for the center atom, i. The angle 0jik between the center atom i and the

45

other two atoms j and k were calculated. The interparticle force is then found by the

derivative of the three-body potential. The radial force and the tangential force are

calculated separately. The components of the force in each dimension i, j, and k are

found by multiplying the previous calculated tangential and radial forces with respect to

the related vectors in the coordinates. The components of force are then summed for

each particle for all interactions.

Updating of the Atom Positions

After the forces on each atom are accumulated from either the two-body or the

three-body potential, the change of positions of particles are updated using Verlet's or

Gear's algorithm at each time step. The details of the Verlet's algorithm59 and Gear's

algorithm60 will be illustrated as follows.

Verlet's algorithm

The advantages of Verlet's algorithm59 are speed and small error. The

derivation for this algorithm is as follows: According to Taylor series expansion, the

position at time (t + At) calculated from time t is

Xi(t+ At)= X(t)+ (dXi(t) )At+ 1 dX(t)) At2 + 1 d3Xi(tAt3 + O(At4);

t 2/ dt2 ) 3! dt3

(3.27)

while the position at (t At) calculated from time t is

Xi(t-At)=Xi(t) d(d1 (t) A) -3At- 1 d3Xi(t)A3+

-' 2 t 3 t3 -+O(At4)"

(3.28)

Adding these two expansions, we can get the new position:

Xi(t +At) =2Xi(t)-Xi(t-At)+( d2x-(t)At2 + O(At4), (3.29)

Where i = 1, 2, 3 denote the Cartesian dimensions of space,

t is a particular time,

Xi(t) is the i Cartesian coordinate at time t,

Xi (t At) is the i Cartesian coordinate at time t At;

At is the length of a time step, and

O(At4) is the remainder terms.

Gear's algorithm

Gear's algorithm60 is another choice to calculate the change of positions of

particles due to two-body interaction. The advantage of Gear's algorithm is that it

gives less energy fluctuation; while the big disadvantage is that it is very time-

consuming. Gear's algorithm predicts molecular positions Xi at time (t + At) using a

fifth-order Taylor series based on positions and their derivatives at time t. Thus, the

derivatives X(i), X(, Xi, Xi(iv), Xi(v) are needed at each step; these are also

predicted at time (t + At) by applying Taylor expansions at time At.

X (t + At) = Xi (t) + X(i) (t)At + Xi(i) (t) (At)2 + Xi() (t) (At)3

2! 3!

+xi(iv) (t) (4t + Xi( (t)(t)5 (3.30)

4! 5!

Xi(i) (t + At) = Xi0) (t) + Xi(ii) (t)At + Xi(iii) (t) (At)2 + x(iv) (t) (At)3

2! 3!

+Xi(v) (t) (At)4 (3.31)

4 !

2 3^

1 tAt+X~iv)() (At X~()(t)(At) (3.2

X(ii) (t + At) = Xi(ii) (t) + Xi(iii) (t)At + Xi(iv) (t) (At + X(v) (t) (3.32)

2! 3

Xi(i") (t + At) = Xi(^) (t) + Xi(iv) (t)At + Xi(v) (t) 2-- (3.33)

Xi(iv) (t + At) = Xi(iv) (t) + Xi(v) (t)At (3.34)

Xi(v) (t + At) Xi(v) (t) (3.35)

In order to correct the predicted positions and their derivatives using the

discrepancy between the predicted acceleration and that given by the evaluated force.

The force at (t + At) obtained from Newton's second law can be used to determine the

acceleration X(ii)(t + At). The difference between the predicted accelerations and

evaluated accelerations is then formed,

AXi(ii) =Xi(ii) (t + At) Xi(ii)P (t + At). (3.36)

In gear's algorithm for second-order differential equations, this difference term

is used to correct all predicted positions and their derivatives. Thus,

Xi = X +ao0AR2 (3.37)

X(i)At = xi(i)PAt + alAR2 (3.38)

X X ") (A t)2

x(ii) (At)2! -x P (At + .X2AR2 (3.39)

2! 2!

xi(iii) (At)3 xi(iii)P (At)3

3!(ii 3!~' + c3AR2 (3.40)

3! 13!

(AtV'(_t_4

xi(iv) (A = X (V)P (At) + a4AR2 (3.41)

4! 4!

xi(V) (At)5 = Xi(V)P (At) + 5AR2 (3.42)

5! 5!cAR

l. A AY()(At)2

where AR2 = AXi' (At) 2

2!

The parameter aci promotes numerical stability of the algorithm. The cXi

depends on the order of the differential equations to be solved and on the order of the

48

Taylor series predictor. Gear determined their values by applying each algorithm to

linear differential equations and analyzed the resulting stability matrices. For a q order

predictor, the values of the ci were chosen to make the local truncation error of

O(Atq+1). Values of the aci for fifth-order predictors are UQ = 3/16, C1 = 251/360,

aC2 = 1, a3 = 11/18, ca4 = 1/6, oa5 = 1/60.

Temperature Calculation

The temperature of the system is controlled by assuming that the system

behaves as a classical statistical system, with average particle kinetic energies given by

the equipartition theorem. The temperature of the system is then calculated by solving

the expression Ek = 3kT/2 for T. The temperature is then adjusted by scaling the

velocities of all particles to give the proper average kinetic energy.

M.D. Video Presentation

Information is collected from the simulation. Any aspect of the simulation can

be recorded, but there are a few which are commonly recorded. These are periodic

snapshots of the structures, instantaneous or average thermodynamic properties, pair

correlation functions, and bond angle distribution functions.

Snapshots of the structure of the system are made by periodically saving all

positions of particles of the system in a data file. These snapshots of the structures

may then later be used to perform a free volume sphere analysis on the structure to

make graphic representations of the structure.

Thermodynamic properties, such as system temperature, energy, and pressure

are calculated. Average of these values may also be kept in a record file. Thus either

instantaneous thermodynamic properties or average properties may be made.

Potential Energy of the System

The potential energy for a system is calculated in a way similar to the

calculation of forces on individual particles. When the neighbors-list is updated,

interaction potentials are found using the potential energy table. The potential energies

for all interactions are then summed to find the potential energy for the system. The

system potential energy is added to the system kinetic energy as previously obtained to

find the system internal energy.

System Pressure Calculation

The system pressure or stress tensor is calculated using the Virial of Clausius,

which has the form

P=p.k.T t-ikT>(j VD(r)) (3.43)

3J J>i

When converted for use in simulations, it takes the form

P=-3 miv +ZiYj Fij, (3.44)

where P is the system pressure in Kbar,

V is the volume of the cell in A3,

v is the velocity of the particle in A/O'10-14s,

F is the interatomic force in 10-12 erg/A,

r is the interatomic separation in A, and

mi is the mass of the particle in 10-24 g.

Pair Correlation Function

Pair correlation function distributions are accumulated when updating neighbors

lists. These are done by incrementing arrays existing for each type of interaction in

which the elements of the arrays correspond to a segment of distance every time a pair

with some separation is found. The contents of the arrays are then converted into

distributions by normalizing using the function

2-nij(r).V

G(r) =- (3.45)

N ij 4rtr2 'Ar

where G(r) is the pair correlation function,

r is the separation distance,

Ar is the separation distance,

V is the volume of the cell,

Nij is the number of pairs found between the distances r and

Ar from the central atom,

Nij m. ni-. nj, (3.46)

where m is the number of accumulations,

ni is the number of atoms of type i,

nj is the number of atoms of type j, and

47rtr-Ar2 is recognized as the volume of the shell in which the non-

central atom may be located about the central atom.

This value is output to a normalized distribution to a data file.

Bond Angle Distribution

Bond angle distributions are found by saving the atom identification numbers of

the coordinating atoms about central atoms during updates of the neighbors lists, then

applying the definition of a dot product to find bond angles. The formula for doing

this is

0 = Cos-1 xxijxik + YijYik + ziJzik (3.47)

9=o--1 k Jl-- y- (3.47)

rij rik

where 0 is the bond angle,

i denotes the apex atom,

j, k denotes the end atoms,

xij, yij, zij are the vector components between i and j,

xik, yik, zik are the vector components between i and k,

rij is the distance between atoms i and j, and

rik is the distance between atoms i and k.

Periodic boundary conditions were used. An energy distribution function is

used to assign the energy (velocity) of each atom at the beginning. The temperature of

the sample can be varied by scaling the velocity of each particle with an appropriate

factor greater than one and letting the total kinetic energy of the particles be equal to

the enthalpy of the desired temperature. The temperature of the sample was kept

constant by scaling the velocity every time step to avoid heating (or cooling) effects.

Introduction of Strain

The determination of the stress-strain behavior will provide an insight into

material behavior for different loading conditions. Uniaxial strains will be applied by

scaling the x-component of all particle positions using the relation

xi=xoi(1 +Li/Sx), (3.48)

where xi and x0i are the scaled and initial x-components of the positions of

particle i, respectively;

Sx is the x-dimension of the cell at the time of scaling, and

Li is the desired incremental expansion factor of the cell in the x-

direction.

Various strain rates can be applied by varying both the value of Li and the

frequency with which the positions are scaled.10 In the calculation, the y and z

directions will be kept constant, the change in volume is directly proportional to the

change in the x direction, Sx.10

CHAPTER 4

EXPERIMENTAL AND SIMULATION PROCEDURES

Experimental Procedure

In this section, the experimental procedure is presented. Sample preparation,

fracture surface analysis, toughness measurements, and fractal dimension determination

will be given in detail.

Sample Preparation

Single crystal silicon was provided by AT&T and IBM. The Laue back

scattering method was used to determine the orientation of the crystal. Several low

index planes were chosen to be the desired fracture planes. They are the {100}, {110}

and {111} planes. After the orientation was determined, single crystal silicon was cut

into flexural bars with the desired orientations. The surface of the specimen was

polished to a l1pm finish.

A controlled flaw was introduced in the desired crystalline plane using a

Vickers diamond pyramid indentation at the center of the tensile surface on flexure

bars. Thus, the flaw was oriented perpendicular to the longitudinal axis of the

specimen. For the {110} orientation, 0.7, 0.9, 1.3, 1.5, 3, 4, and 5 Kgw indentation

loads were applied to produce a controlled flaw in each bar. For the {100} and {111}

orientations, 1 and 2 Kgw indentation loads were applied. Figure 4.1 shows a

schematic arrangement to demonstrate the location of the indentation on the tensile

surface. Flexural testing was performed on an Instron testing machine in air at room

temperature. Loads to failure were recorded. Fracture surfaces were examined to

(a) {1 00} fracture plane

100}

(b) { 11 O} fracture plane

{110}

{11}

{111} fracture plane

{ 10}(

{I12} '{111}

S {112}

110o} {111}

Figure 4.1 Three orientations are chosen for comparison.

insure the samples had failed at the indent site. The stress to failure was calculated

from the dimension of each sample:

a 3PL (4.1)

where Ga is the fracture strength,

P is the load to failure,

L is the load span,

b is the width of the bar, and

h is the height (thickness) of the bar.

Fracture Surface Analysis

The procedures for fracture surface analysis have been well described in the

literature.6 The size of the critical flaw was determined by modeling the flaw as an

idealized elliptical crack of depth a and half width b. The equivalent semi-circular

crack radius, c, is determined by c = ab and is shown in Fig. 4.2.30

The fracture surfaces were observed using a light microscope to locate the

flaws. Failure origins were located by observing particular fracture markings on the

fracture surfaces.61 Crack branching and mirror features were observed and were

used to locate the flaw.

Toughness Measurements

From the measurement of the critical flaw size, c, and fracture strength, aa, the

critical stress intensity factor, K., for those with residual stress caused from indentation

can be determined using

Kc = 1.65ac12. (4.2)

(b)

/a/

2b

Fig. 4.2 Micro-indentation cracked bar under bending. (a) A schematic diagram

of Vickers indentation on tensile surface. (b) An elliptical flaw on the fracture cross

section.

Furthermore, we can compare the obtained results from fracture surface analysis with

those from the strength-indentation technique by using:

Kc = 0.59(E / H)8(CaP13)34. (4.3)

Here, the fracture toughness is given in terms of the critical intensity factor,

Kc, instead of Kic. Kic is the critical stress intensity factor for mode I. Kic is usually

determined from a prescribed fracture mechanics test with no residual stress associated

with the crack. Kc is the resistance to crack propagation in the presence of a local

residual stress. Kic and Kc are expected to be close in value. During the test, the

mode one condition is not guaranteed because the fracture paths tend to go to the

easiest cleavage plane. Also, the stress condition on the fracture surface is not in a full

tension mode. Indeed, sometimes it is a combination of mode I and mode II. Thus, Kc

is an approximation of the mode I fracture toughness of the material.

Fractal Dimension Determination

Single crystal silicon fracture samples were carefully cleaned and coated with

nickel. The samples were then potted in epoxy. In general, there are two ways to get

contour lines: one is to polish the sample approximately parallel to the fracture surface

and the other is to polish the sample approximately vertical to the fracture surface.

Thus, the so-called horizontal contour (fractal) dimension or vertical profile dimension

can be obtained. The samples in this study were polished approximately parallel to the

fracture surface, as shown in Fig. 4.3. As the fracture surface is first encountered, a

section of the fracture surface appears in the polishing plane. These sections appear as

islands in the polishing plane, and are called slit-islands. The nickel coating performs

two functions: It provides good contrast during polishing and it helps to hold the

fracture surface together during polishing. As polishing proceeds, these islands begin

potted in epoxy and

polished parallel to surface

top view

Fig. 4.3 Fractured samples are encapsulated in epoxy and polished approximately

parallel to the fracture plane. (Top View) Islands emerge in the polished plane.

fractured bar

to grow. The perimeter of the islands presents a line or section of a line that can be

measured according to Richardson's equation:

L(S) = k S1-D, (4.4)

where L(S) is the length of a section of a line along the perimeter of the slit

island and its value depends on the measuring scale, S;

S is the measuring scale used to measure the perimeter and its value

ranged from 5 to 100 utm;

D is the horizontal contour (fractal) dimension; and

k is a constant.

Richardson plots detail the change in measured length of a line as a function of

scale. Construction of a Richardson plot for a fracture surface requires a line that is

representative of the fracture surface and a range of scales for measuring that line.

At the first emergence of an island, polishing must be carefully performed. The

surface was polished to a 1 gm finish. Since the epoxy is transparent, the exact

location on the fracture surface can be observed, so that measurements can be taken for

selected segments of the perimeter, e.g., in the crack branching region. Polaroid

photographs were taken at a magnification of 400x and combined in a montage. This

montage was then measured with dividers set to various openings, as depicted in Fig.

4.4. The length of the section of perimeter was measured for each divider setting. In

this way, the line or perimeter length was computed as a function of scales. A log-log

plot of length vs. scale results in a line with slope equal to 1-D.

The scale is a measure of discernibility. As the scale becomes finer and finer,

we observe greater and greater detail. If the scale is larger than the largest features,

the observation is insensitive to those features. A curve will begin to look fractal only

after the scale becomes smaller than such features. The scale lengths used in this study

ranged from 1 unit length to 16 unit lengths (1 unit length is about 5 pin).

A montage of the perimeter is measured with dividers.

Fig. 4.4

M.D. Simulation Procedure

Before performing the MD simulation, parameters24 for the potential and initial

conditions should be given beforehand. Two data files which contain these data are

needed in order to perform the MD simulation. One is called the input file, which

includes the initial conditions and the parameters of the potentials used. The other one

is the atom-position file which contains the initial position of each Si atom.

Several tasks were set to be accomplished using MD simulation. One was to

compare the toughness of Si for different orientations. The difference of stress-strain

curves for different strain rates could be obtained to compare the toughness as a

function of strain rate. Another was to obtain the fractal dimensions for different

fracture planes in order to compare them with the experimental results. Another goal

was to obtain the fracture strengths for each simulation. The comparison of fracture

strength due to different strain rates has been performed. The comparison of fracture

strength due to different crack sizes has also been performed.

Determination of Input Data

Before performing the MD simulation, several conditions should be decided

first. All those conditions should be illustrated in the input file for the MD simulation.

Those conditions are the desired orientation, temperature, strain rate, length of each

time step, length of the simulation, and the adiabatic or isothermal state.

At first, after choosing the desired orientation, a data file which contains the

number and the positions of those atoms should be constructed. The data file will

determine the initial position of each atom for the MD simulation. In this study, the

atom position for the ideal crystal structure will be given at first and the crystal will be

allowed to thermally vibrate for one unit time, 1 pica-second, before performing strain

pulling.

The temperature to perform the test should also be given. The chosen

temperature is important because it determines the energy state of the system. Here the

temperature or initial temperature is chosen to be 300 K.

The strain-rate should be given, also. Different strain-rates will give different

results. Basically, two categories will be given. One is higher than the speed of sound

while the other is slower than the speed of sound. Generally, results from different

strain-rate experiments will be compared for the same orientation. While for different

orientations, a strain rate of 0.2 was usually chosen. The reason for selecting 0.2 will

be explained in the results and discussion section.

A careful choice of the length of each time step is very important. It plays a

very important role not only on the stability of the system but also on the length of the

simulation time. If the length of the time step is chosen to be too long, atoms will

move too far in each time step. Thus the simulation will either be less realistic or

become unstable. If the time step is too small, atoms just move a very small distance at

each time step and waste a lot of computation time. For silicon, 0.005 pica seconds is

a good choice for the Stillinger-Weber and (ionic) Coulombic potentials because the

thermal vibration period is found to be 0.0766 ps for Si.24 The length of the

simulation is determined by the fracture of the specimen. For a lower strain-rate test, a

longer time should be given.

The choice of an adiabatic or isothermal condition is very important to do MD

simulations. If the adiabatic condition is used, the temperature will increase several

thousands Kelvin degrees and the system becomes unstable after fracture. If only the

Stillinger-Weber potential is used, the temperature after fracture is about 7,000 K. If

the ionic Coulombic potential is added in, the temperature after fracture will be about

50,000K. The reason is that the ionic Coulombic potential is always repulsive. Thus,

the energy state of the system will increase as more positive charged Si atoms are

produced during the fracture process. On the other hand, the isothermal condition will

result in a more stable thermal state for this study since the temperature remains

constant. Thus, the isothermal condition was used in this study. The temperature vs.

strain curves for both isothermal and adiabatic condition are shown in Fig. 4.5 and 4.6.

The combinations of the previous conditions generally should satisfy different

needs. After these conditions have been decided, one can begin to perform the MD

simulation.

Simulation Procedure

At first, toughnesses for different orientations were compared. Three different

orientations were chosen in this study. The fracture planes were chosen to be { 100},

{110} and {111} as shown in Fig. 4.7. A Periodic boundary condition was applied

during the simulation. A constant expansion rate in the x-direction was applied for

each simulation.

The stress-strain curves for different strain rates were compared in order to find

the dependence of strength on strain rate. The strain rates were chosen to be 0.1, 0.2,

0.5, 1, 2, 5 x-length/ps.

In the Griffith criterion the strength of brittle materials depends on the crack

size at which fracture begins. The introduction of a crack in the ideal crystal will cause

a decrease in the fracture strength. Larger cracks are expected to result in lower

strengths. A crack can be introduced in two ways. One type of crack simulates a void

by removing a cluster of atoms. The other type simulates a planar sharp crack by

removing a layer of atoms.

315

.. 310

305

S300

g 295

290 I

0 0.2 0.4 0.6 0.8 1

strain

Fig. 4.5 The temperature vs. strain curve for isothermal condition. The

temperature remains at about 300K before the moment when fracture takes place at

strain 0.45. After then the temperature goes back to about 300 K.

10000

9000

.' 8000

^ 7000 -

S6000 / ^

S5000

S4000

S3000

S2000

1000

0

0 -- ---- )--------

0 0.2 0.4 0.6 0.8 1

Strain

Fig. 4.6 The temperature vs. strain curve for adiabatic condition. The

temperature remains at about 300K before the moment when fracture takes place at

strain ; 0.45. After fracture, the temperature of the system rise very fast up to 7000

K, then slows down.

fil

Fig. 4.7 The 11001, 11O0, and {11} planes are those fracture planes chosen in

the MD simulation.

/ :/ /no //{i

/{^\ /_^- /J -- J

Fig. 4.7 The {100}, {110}, and {111} planes are those fracture planes chosen in

the MD simulation.

Fractal Analysis Using Simulation Results

After the fracture has occurred in the simulation as shown in Fig. 4.8, an

important objective is to obtain the generated fracture surface from the MD

simulation.Each atom is assumed to be a sphere with an electron cloud of a prescribed

radius. The radius of each atom is assumed to be the radius of the potential field.

Thus the fracture surface appears as a surface with intersecting spheres, as shown in

Fig. 4.9. After obtaining the surface, a contour plane is used to intersect the spheres to

obtain a "slit-island." The obtained slit-island appears as a plane with lots of circular

disks on it as shown in Fig. 4.10. The slit-island technique as described in the fractal

dimension determination section in the experimental procedure was applied here to

obtain the fractal dimension. The perimeter of the selected island was measured using

different scales as described before and the fractal dimension can be obtained from the

Richardson plot.

The stereo-pair picture shows a frozen frame during fracture.

Fig. 4.8

Fig. 4.9 If each atom is assumed to be as a potential field with sphere shape. The

fracture surface will look like a surface with intersecting spheres.

.. . ..

~ __ ........ .,

I4e

Fig. 4.10 The slit-island obtained from the simulated fracture surface looks like a

plane with lots of circular disks on it.

CHAPTER 5

RESULTS AND DISCUSSION

Fracture toughnesses and fractal dimensions were obtained in two ways, one

was from experimental measurement and another was from molecular dynamics (MD)

simulation. Fracture toughnesses and fractal dimensions were investigated for three

different fracture planes. This chapter discusses the toughness results, the fractal

dimension results, and how the fractal dimensions for each orientation are related to the

toughness results. This chapter is divided into three sections: experimental results,

simulation results and comparison between measured and simulation results.

Relationships between variables will be discussed in each section.

Experimental Results

In this section, results obtained from the experimental measurements will be

given. These include fracture surface analysis, fracture toughness measurements,

fractal dimension measurements, and the relationship between fractal dimension and

fracture toughness.

Fracture Surface Analysis

After the specimens were broken, the fracture surfaces were examined. The

critical flaw size was measured to calculate the critical stress intensity factor using Eq.

(4.2). As in Fig. 4.1, three fracture planes, i.e., {100}, {110} and {111}, were chosen

and two different tensile surfaces in the {110} and {111} fracture planes were tested.

Thus, there were a total of 5 different groups of data. The fracture surfaces are shown

in Figs. 5.1, 5.2, 5.3. As shown in Fig. 5.1, the fracture surface of the {100} fracture

plane appears to be the most tortuous one. It has the smallest mirror region when

compared with the other two orientations. The fracture surface of the {110} fracture

plane with a {100} tensile surface has a "Batman"-shaped mirror region as shown in

Fig. 5.2(a). The fracture surface of the {110} fracture plane with a {110} tensile

surface has an inverted volcano-shaped mirror region as shown in Fig. 5.2(b), while

the fracture surface of the {111} fracture plane as shown in Fig. 5.3 is relatively

smooth compared with the other orientations. "River-marks", i.e., twist hackle or

cleavage marks, appear on most of the fracture surfaces of this fracture plane. "River-

marks" are twist hackle fracture features which appear to spread apart as the crack path

moves away from the origin of fracture, thus "pointing" back to the origin. The {111}

fracture plane is the easiest cleavage plane, as reported before.22 Shown in Fig. 5.3(a)

is the fracture surface of the {111} fracture plane with a {110} tensile surface while

that in Fig. 5.3(b) is the fracture surface of the same fracture plane but with a {112}

tensile surface. More data are listed in Tables 5.1, 5.2, and 5.3.

For fracture on the same plane, there should be no structural difference during

the fracture process. However, the loading, i.e., strain would be expected to be

different in different locations due to elastic anisotropy and a changing strain field,

which in turn is due to the type of strain applied here, i.e., bending. The fracture

surface can appear different due to a change in the tensile surface in the bending test.

Elastic anisotropy is the reason the fracture surfaces look different in the same fracture

plane without a structural change. Thus the fracture surfaces of the {110} fracture

planes as shown in Figs. 5.2(a) and (b) look different because the tensile surfaces are

different. An interesting result occurs if we rotate the fracture surface by 90 degrees

with a {110} tensile surface, as shown in Fig. 5.4; a half mirror boundary of this

fracture surface is similar to a half mirror boundary of the fracture surface on the

A typical fracture surface on the {100} fracture plane.

Fig. 5.1

{100} }

(b)

MOW,

_- d.. JB &. .

{110}

^--{

Fig. 5.2 (a) A typical fracture surface of the {100} tensile surface on the {110}

fracture plane. A "Batman"-like mirror is obvious on the fracture surface. And a

typical flaw is easily seen on the fracture surface. (b) A typical fracture surface of the

{110} tensile surface on the {110} fracture plane. An inverse volcano-like mirror is

obvious on the fracture surface, and a typical flaw is easily seen on the fracture

surface.

W

(a)

{11} }

(b)

F{11 2}

Fig. 5.3 (a) A typical fracture surface of the {110} tensile surface on the {111}

fracture plane. (b) A typical fracture surface of the {112} tensile surface on the {111}

fracture plane. It is found that both fracture surfaces look similar even though the

tensile surfaces are different.

Table 5.1 Data for Si samples fractured in the {100} plane

{100} fracture plane ,<110> tensile surface _______________

indent flaw load KIc* KIc**

S load widththickness b a size P strength 0.59(E/ 1.65*

No. Kg mm mm mm mm mm lb MPa MPam^l/2 MPam^l/2

al 0.5 5.70 5.35________

a2 0.5 6.75 5.00 0.068 0.068 0.068 96 102.5 1.31 1.40

a3 0.5 5.85 4.95 0.080 0.080 0.080 65 81.7 1.10 1.21

a4 0.5 6.70 5.40 0.056 0.052 0.054 118 108.8 1.37 1.23

_____average 1.26 1.28

--___ -- dev 0.11 0.08

bl 1 6.25 5.05 0.127 0.094 0.110 63 71.2 1.18 1.23

b2 1 5.90 5.25 0.080 0.080 0.080 78 86.4 1.37 1.28

b3 1 5.80 5.05 0.124 0.124 0.124 54 65.8 1.11 1.21

b4 1 5.90 5.00 0.114 0.072 0.091 63 77.0 1.25 1.21

b5 1 5.80 5.10 0.127 0.099 0.112 60 71.7 1.19 1.25

_____average 1.23 1.24

_____ dev 0.09 0.03

cl 2 .6.25 5.20 0.201 0.132 0.163 55 58.6 1.22 1.23

c2 2 5.90 5.15 0.159 0.131 0.144 57 65.6 1.32 1.30

____. average 1.27 1.27

___________dev 0.05 0.03

t____otal average 1.24 1.26

____________total dev 0.09 0.06

__three-point bending span=27mm_________________

KIc is calculated from Eq. (4.3) **KIIc is calculated from Eq. (4.2)

Table 5.2 Data for Si samples fractured in the {110} plane

{100} tensile surface ___ _____________________

indent flaw load ____ Kc* K 1c**

No. ToadT width thickness b a size P strength 0.59(E/ 1.65*

SKg mm mm um um urn N MPa MPam^l/2 MPamAl/2

1 0.9 7.44 7.03 66.2 84.8 74.0 322.6 103.1 1.46 1.46

2 0.9 7.31 7.02 70.2 66.3 68.0 0.0. 107.4 1.51T 1.46

3 0.9 7.31 7.00 U 66.2 79.5 72.0 249.2 81.8 1.23 1.14

4 0.9 73T 6.99 79.5 92.8 86.0 262.5 86.7 1.29 1.32

-53 0.9 7.53 7.15 -95T 98.0 96.5 266.9T 81.5 1.23 1.31

6 0.9 6.87 7.16 92.8 90.0 91.5 238.0 79.4 1.21 1.25

7 0.9 6.58 7.03 94.0 94.0 195.8' 70.8 1.10 1.13

8 0.9 6.89 7.01 97.0 97.0 242.6 84.2 1.26 1.35

9 0.9 7.07 7.07 102.0 f102.0 202.4 67.3 1.06 1.12

10 0.9 6.81 7.39 953. 89.8 89.8 233.6 73.9 1.14 1.15

11 0.9 75735 6.95 89.8 92.4 198.0 63.9 1.02 1.01

12 1.4 7.11 7.37 130.0 143.0 136.3 164.6 50.1 0.95 0.95

13 1.4 6.88 7.11 136.0 80.87 __

14 1.4 7.84 7.16 32.5 124.0 127.9 193.5 56.6 1.05 1.05

5 1.4 T 7.07 7.44 119.0 106.0 112.1 229.1 68.8 1.21 1.23

16 1.8 6.83 6.92 126.0 172.0 147.0 198.0 71.1 1.32 1.41

17 1.8 7.43 02242.5 77.8 1.41 _____

18 1.8 7.82 6.94 220.0 ___ 220.2 68.7 1.29 _____

19 1.8 6.9 7.35 185.0 185.0 200.2 63.2 1.21 1.41

20 0.53 6.87 7.28 63.4 52.8 57.8 269.2 87.0 1.12 1.09

21 0.5 7.36 6.60 352T 52.8 52.8 262-.5 96.3 1.20 1.15

22 2 7.57 6.90 132.0 142.6 137.2 204.3 66.8 1.29 1.29

23 2 7.70 7.32 190.0 158.4 173.5 191.3 543.5 1.11 1.18

24 3 7.41 7.24 T195.4 142.6 166.9 209.1 63.3 1.37 1.35

25 3 6.87 6.95 169.0 184.8 176.7 186.9 66.2 1.42 1.45

26 5 6.29 7.15 322.1 294.6 307.6 133.5 48.8 1.28 1.41

27 5 7.45 7.37 316.8 253.4 283.3 169.1 49.1 1.29 1.36

________ avg. 1.19 1.23~

________________dev. 0.08 0.08

{110} tensile surface _____ _____ ______

indent flaw load Kc1* K~ Klc**

No. load width thickness b a size P strength 0.59(E/ 1.65

- Kg mm mm umr um um N MPa MPam^ 1/2 MPamA1/2

41 0.9 7.35 6.97 T1584 105.6 129.3 160.2 52.76 0.89 0.98

42 0.9 7.55 7.39 95.0 95.0 95.0 278.1 79.30 1.20 1.26

43 0.9 7.10 7.56 89.8 100.3 94.9 233.6 67.71 1.07 1.07

44 0.9 6.93 7.84 1056 105.6 105.6 229.1 63.24 1.01 1.06

45 0.9 6.76 7.03 84.5 79.2 81.8 204.7 72.06 1.12 1.08

46 0.9 7.44 7.26 843. 95.0 89.6 215.8 64.73 1.04 1.01

47 5.0 7.08' 7.21 380.1 337.9 '3584 100.1 31.99 0.93 0.99

48 5.0TT739 7.26 269.3 316.8 292.1 142.4 43.00 1.17 1.21

______avg. 1.05 1.07

____ ____dev. 0.10 0.09

Klc is calculated from the Eq. (4.3) ** Klc is calculated from the Eq.(4.2)

Table 5.3 Data for Si samples fractured in the {111 } plane

110 ,tensile surface__________________

indent _____flaw load ____KIc* KIc**

load width thickness b a size P strength 0.59(E/ 1.65*

No. Kg mm mm nmm mm mm lb MPa MPamA^l/2MPam^l/2

cl 1 6.15 4.00 0.082 0.092 0.087 33 76.1 1.20 1.17

c2 1 5.80 4.95 0.093 0.102 0.098 46 73.4 1.17 1.20

c3 1 6.30 5.05 0.091 0.103 0.097 40 56.5 0.96 0.92

A1 1 6.15 5.60 0.110 0.112 0.111 82 70.9 1.14 1.23

A2 1 6.20 5.20 0.103 0.124 0.113 78 77.6 1.22 1.36

_____avg 1.14 1.18

________ dev 0.09 0.15

dl 0.7 6.15 4.65 0.078 0.074 0.076 44 75.1 1.09 1.08

d2 0.7 6.20 5.15 0.058 0.063 0.061 72 99.4 1.34 1.28

d3 0.7 6.10 4.90 0.080 0.086 0.083 48 74.4 1.08 1.12

d4 0.7 6.25 5.05 0.065 0.085 0.075 56 79.7 1.14 1.14

avg 1.16 1:15

_______________dev 0.11 0.07

el 2 6.30 4.80 0.171 0.200 0.185 34 53.2 1.09 1.19

e2 2 6.25 4.75 0.149 0.155 0.152 40 64.4 1.26 1.31

e3 2 5.70 4.70 0.172 0.173 0.172 29 52.3 1.08 1.13

e4 2 5.95 4.95 0.160 0.118 0.137 38 59.1 1.18 1.14

_________avg 1.15 1.19

_________dev 0.07 0.07

{112} tensile surface___________________

fl 2 6.15 6.25 0.108 0.115 0.111 68 64.2 1.26 1.12

f2 2 5.35 6.20 0.125 0.144 0.135 60 66.2 1.29 1.27

13 2 5.25 6.15 0.165 0.125 0.143 50 57.1 1.15 1.13

f4 2 5.80 6.10 0.165 0.183 0.174 48 50.5 1.05 1.10

avg 1.19 1.15

dev 0.09 0.07

gl 1 6.40 5.80 0.079 0.094 0.086 67 70.6 1.14 1.08

g2 1 5.50 6.20 0.105 0.105 0.105 70 75.1 1.19 1.27

g3 1 5.85 5.65 0.086 0.099 0.092 52 63.2 1.05 1.00

B1 1 6.30 6.05 0.109 0.098 0.103 101 73.1 1.17 1.23

B2 1 5.80 5.70 0.086 0.071 0.078 93 82.3 1.28 1.20

avg 1.16 1.16

____dev 0.07 0.10

KIc is calculated from Eq. (4.3) **KIc is calculated from Eq. (4.2)

A,B span = 25.5 mm J c,d,e,f, and span =34

three-point bending I Jf I | n

(a)

(b)

4;

{100} 1 }

/{m /{

-i

K

-~ -~ -

Fig. 5.4 An interesting result is that if we rotate the fracture surface of the one

with a { 110} tensile surface 90 degrees, as shown in (a), then a half mirror boundary of

this fracture surface is similar to half a mirror boundary of the fracture surface on the

{110} fracture plane with a {100} tensile surface without rotation, as shown in (b).

This implies that a fracture mirror boundary for a specimen fractured in tension from

an internal defect would mimic the trace of the elastic constant values in that plane.

{110} fracture plane with a {100} tensile surface without rotation.50 This implies that

a fracture mirror boundary for a specimen fractured in tension from an internal defect

would mimic the trace of the elastic constant values in that plane.52

Freiman et al.52 found that the fracture mirror shapes in sapphire for pseudo-

cleavage type fracture have much more complex geometry, as shown Fig. 5.5, than the

similar fractures in glasses which yield more-or-less circular fracture mirrors. They

relate the mirror boundaries to the critical stress intensity factor. The complex mirror

shapes investigated in sapphire implies that the critical stress intensity factor along the

different crystallographic directions of the mirror boundary is not constant. Young's

modulus also plays an important role for the complex shape of mirror boundary. The

Young's modulus normal to the crack plane is a constant. However, the Young's

moduli either parallel or perpendicular to the mirror boundary in the crack plane can

vary as one traverses around the boundary.62 As seen in Fig. 5.5, The Young's

moduli for directions parallel and perpendicular to the crack front/mirror boundary

plotted as a function of angle of rotation about the tensile axis possesses maxima-

minima patterns which are similar to the fracture mirror geometry. A similar result for

single crystal silicon is expected and a schematic of the expected mirror boundary is

plotted in Fig. 5.5.

In contrast to the fracture surface of the {110} fracture plane, the fracture

surfaces of the {111} fracture plane look the same even when the chosen tensile

surfaces are different as shown in Fig. 5.3. Most probably, the fracture surface of the

easiest cleavage plane will not be affected by the chosen tensile surface because the

crack growth will remain on the easiest cleavage fracture plane as the lowest energy

path. Also, the size of the specimen is not large enough to show crack branching at

these tested fracture stress levels. Thus, the fractography would appear the same no

matter what the change in the tensile surface at these test specimen sizes. The available

fracture paths for the specific fracture planes are another important factor to determine

<1014>

/I -- .

Fig. 5.5 (a) Anisotropy mirror shapes on fracture surfaces in single crystal

alumina; the corresponding relative magnitudes of Young's modulus for directions (in

the plane) parallel and perpendicular to the local mirror boundary are given. Larger

values of E are reflected as greater distances from the origin. (b) A schematic of the

expected mirror boundary for the {110} fracture plane is plotted.

the topography of the fracture surface. A more detailed discussion of macroscopic

crack branching follows.

The crack branching plane is an interesting topic in this study. From the

observation of the obtained fracture surfaces, the fracture surface of the { 111 } fracture

plane is relatively smooth such that no crack branching is observed as shown in Fig.

5.6. For fracture in the {110} plane, the crack either doesn't branch or tends to branch

to the available {111} planes as shown in Figs. 5.7 and 5.8. In Fig. 5.7, the crack

growth (shown by shaded regions) occurs on the {110} plane first and deviates to

{111} plane as shown in (c) and (d) of Fig. 5.7. In Fig. 5.8, the crack growth (shown

by shaded regions) occurs on the {110} plane first and deviates to {111} planes as

shown in (c) and (d) of Fig. 5.8.

Crack branching for the {100} fracture plane is complex. As seen in Figs. 5.9

and 5.10, the system of crack branching angles varies for one fracture plane as well as

for the loading direction changes, e.g. from the {100} tensile surface to the {110}

tensile surface.

For those fracture surfaces of the {100} fracture plane with a {100} tensile

surface, the crack easily branches to the {111} plane as shown in Fig. 5.9. If the

tensile surface is chosen to be {110}, the crack will remain on the {100} plane for

about 2 to 3 mm before branching to the {111 } plane if the indentation load is larger

than 2 Kgw, as shown in Fig. 5.10(b). If the indentation load is smaller than 2 Kgw,

the crack will branch in a distance smaller than 2mm as shown in Fig. 5.10(c).

It is also found that the fracture mirror region is not really very smooth.12 The

reason that we say it is "mirror-like" is because the tortuosity of the region is smaller

than the wave-length of visible light as shown in Fig. 5.11. Indeed, if a scanning

tunneling microscope (STM) is used to observe the mirror region, then the surface is

not smooth at the scale at which the STM is used. 12 As in Fig. 5.12, the surface on

the mirror region is quite rough at that magnification. Thus, the geometry in the

,, {110}

{io

,{111}'{

, ;11 } L

i{\ {11:

{111} {111}

I :i!l 0!

Crack path for the { 111 } fracture plane.

Fig. 5.6

{low,

1110}

fioo

Fig. 5.7 The crack path for the {110} fracture plane with a {100} tensile surface.

The crack growth (shown by shaded regions) occurred on the {110} plane first and

deviated to the {111} plane as shown in (c) and (d). (a) shows a schematic of the

fracture bar. (b), (c), and (d) show the different modes observed.

{111}

{11

! 110

{100} .

I, I10:

1

f110}

*10I 10

*lO01. \ 'r ;

_ [110o

Fig. 5.8 The crack path for the {110} fracture plane with a {110} tensile surface.

The crack growth (shown by shaded regions) occurred on the {110} plane first and

deviated to the {111} planes as shown in (c) and (d). (a) shows a schematic of the

fracture bar. (b), (c), and (d) show the different modes observed.

SI 10k

:110:

>* ,

7 L {100}

10

{l00}

Fig. 5.9 The crack paths (shown by shaded regions) for the {100} fracture plane

with a {100} tensile surface. (a) shows a schematic of fracture bar. (b) shows

different modes observed.

S {11O}

101

S{110}

Fig. 5.10 The crack paths (shown by shaded regions) for the {100} fracture plane

with a {110} tensile surface. (a) shows a schematic of fracture bar. If the tensile

surface is chosen to be {110}, the crack will remain on the {100} plane for about 2 to 3

mm before branching to the { 111} plane if the indentation load is larger than 2 Kgw as

shown in (b). If the indentation load is smaller than 2 Kgw, the crack will branch in a

distance smaller than 2 mm generally as shown in (c).

mirror 4-

region

mist

region

-M _- hackle

region

wavelength

of light

_4

Fig. 5.11 Schematic of vertical profile of the mirror, mist and hackle regions of a

fracture surface. The reason why the mirror region is mirror-like is because the tortuosity

of the region is smaller than the wave-length of visible light as shown above.

Fig. 5.12 The fracture surface obtained from scanning tunneling microscopy.

(Courtesy of Mitchell and Bonnell12)

mirror region appears to be smooth to the unaided eye, but appears to be tortuous if the

STM is used as the tool of observation.

The reason that the fractographic features described above are important is that

the distances to crack branching have been shown to be related to different energy

levels at branching and the fracture toughness of the material.50 In addition, fracture

surface features have also been shown related to the fractal dimension of the surface.63

Thus, the fracture surface features are intimately connected to the fracture process as

described by the fracture toughness, branching constants and the fractal dimensions.

These inter-relationships will be discussed later in this chapter.

Fracture Toughness Measurement

Single crystal Si is anisotropic and the fracture toughness varies with the

fracture plane (and the tensile surface if a bending technique is used to fracture the Si).

Chen and Leipold22 determined the toughness for single crystal silicon tested in

different orientations. They found that the {111} fracture plane is the easiest cleavage

plane and the {100} fracture plane is the toughest one. They used an indentation

technique to control the crack sizes on the tested specimens before fracture. However,

in their work, they did not consider the residual stress caused by the indentation. The

equation which they used to calculate the critical stress intensity doesn't consider

residual stress:

Kc = YMB(na/Q)l/2, (5.1)

where Kc is the critical stress intensity factor,

cy is the maximum tensile stress,

a is the flaw size,

MB is the elastic stress intensity magnification factor, and

Q is the flaw shape parameter calculated from the elliptical integral of

the second kind.

From the given flaw size data in their paper,22 the MB value was calculated to

be 1.03 and Q was approximately 1.47. After substituting these data into equation

(5.1), the equation can be approximately rewritten as

Kc = 1.24oc/. (5.2)

If we modify Eq. (5.2) to include the effect of residual stress (due to the indentation

1/2404

process), the new modified equation becomes Kc = 1.65a fc i.e., Eq. (4.2).4044

Their data using Eq. (4.2) fall in the same range as those measured in this study. The

measured data and modified data from Chen and Leipold are presented in Table 5.4.

Table 5.4 Comparison of toughness values using different techniques

S.I. F.S.A. literature

fracture plane tensile surface toughness* toughness** toughness#

______ __________(MPaim) (MParim) (MPaJfm)

{100} {110} 1.24 0.09 1.26 0.06 1.26 0.07

{110} {100} 1.19 0.08 1.23 0.08 1.19 0.13

{110} 1.05 0.10 1.07 0.09

{111} {110} 1.19 0.10 1.17 0.09 1.09 0.09

{112} 1.21 0.09 1.16 0.08

*Calculated from the strength indentation (S.I.) technique,

KC =0.59E 1/8(afpl1/3 )3/4

1/82

Kc= 0.59 (fP1/3 .

**Calculated from fracture surface analysis (F.S.A.), Kc = 1.65afc12.

#By Chen and Leipold.22 Originally these values were calculated without considering

residual stress and were found to be 0.950.05 ({100} plane), 0.900.11 ({110}

plane), and 0.82+0.07 ({111} plane) MPaVm.

Gilman2, Jaccodine,64 Myers & Hillsberry,65 St. John,57 and Kalwani66 also

investigated the fracture toughness in single crystal silicon. Those data and applied

techniques are listed in Appendix A.

Single crystal Si is a diamond cubic structure which is considered to be one type

of face centered cubic structure. A simple approach67 which has been accepted to

estimate the surface energy of different planes in a FCC structure is to calculate the

number of dangling bonds on the specific plane for the FCC materials. The surface

energy can be obtained using the equation:

y = N "Nhkl} n{hkl}, (5.3)

2

where e is the energy needed to fracture the Si-Si bond,

is used because two surfaces are generated,

2

N{hkl} is the number of atoms per area on the {hkl} plane, and

n{hkl} is the number of dangling bonds per atom on the {hkl} plane.

By using the given equation, we can find the surface energy for the FCC

structure on the {100}, {110}, and {111} planes separately as shown in Fig. 5.13.

Since the number of atoms on the {100} plane is 1 + 1/4 4 as shown in Plot (b),

which is 2, and the area A where those atoms lie is a a where a is the length of the

side of a FCC unit cell, N{100} is found to be 2/a2. Since the number of dangling

bonds per atom on the {100} plane is 4 as shown in Plot (c), Y{100oo} is found to be

1 2 1 2

2 -E .- 4. In the same way, Y{1o10} is found to be 2- E 2 5, and Y{11 is found to

2 a-;2 x/2a-

1 2

be .6. .-3. Thus y{100} : Y{ll0} : Y{lll} = 8 : 7.07 : 6.92 = 1.16 : 1.02 :1.

2 "r3a 2/2

This ratio means that the {111} plane for the FCC structure has the lowest surface

energy and the {100} plane has the highest surface energy.

(100) = a

A(1O) 2a2

~A(=32

A(llIl)=Vf3a /2

#. of atoms

= 4x1/4+1

#. of atoms

=4x1/4+2x1/2

#. of atoms

= 3x1/6+3x1/2

Fig. 5.13 The surface energy for FCC structures on the (100), (110), and (111)

planes can be found using a dangling bond calculation.

n(loo)= 4

00

00

n(110)= 5

n(11)= 3

However, single crystal Si has a diamond cubic (DC) structure which, in fact, is

different from FCC because it is not closed packed. We can find the surface energy for

the DC structure on the {100}, {110}, and {111} planes separately as shown in Fig.

5.14. However, the DC structure is not a close-packed structure, and the dangling

bonds will not be linked to all neighbors. Only some specific neighbors in the first

shell will share the dangling bond. The previous method cannot be completely

followed and some modifications must be made. As shown in the figure, the dangling

bond is expressed as a line. The surface energy can be obtained using the equation:

1 1 (5.3)

Y=-2"{hkl} "A' A5'3)

2A

where 6 is the energy needed to fracture the Si-Si bond,

1.

is used because two surfaces are generated,

2

n{hkl} is the number of dangling bonds on the {hkl} plane, and

A is the area where those atoms lie on.

Since the number of dangling bonds on the {100} plane is 4 as shown in Plot

(c), and the area A where those atoms lie is a a and a is the length of the side of a

silicon DC unit cell. y{100oo} is found to be I .E If s = 3.5-10-12 erg/pair-atom22

2 a2

and a = 5.43 A is substituted into the equation, 7{100oo} is found to be 2.3 J/m2. (The

Si-Si bond energy is found to be varied from 3.2 to 3.6 erg/pair-atom.24'64'68'69) In

1 4 _

the same way, y{110o} is found to be -*1 = 1.67 J/m2, and y{11} is found to be

2 72:a2

1 2

-.e. 1.36 J/m2. These values are close to the experimental value, 2.1

2 73a2/2

J/m2, which is obtained using double cantilever beam tests.70'71 Thus Y{100} : Y{110o}

Y{111} = 2.3 J/m2 : 1.67 J/m2 : 1.36 J/m2 = 1.69 : 1.23 : 1. So, from dangling bond