MD simulation of crack propagation in
brittle crystals 118094  Introduction to Computational Physics  Fall 2012/13 Liron BenBashat Bergman 
Overview  Theory  Setup  Results  Manual 
Simulation Setup
Molecular Dynamics simulation Molecular dynamics (MD) is a computer simulation
technique which enables the modeling of material behavior as a function of time. The physical
principle underlying MD simulation is the BornOppenheimer adiabatic
approximation. Considering a Newtonian behavior, the equation of motion
of a system consisting N particles which apply on each other a force
field,
, or equivalently potential field U is given by:
,
, and
are the
position, acceleration, and the mass of the ith atom respectively and
is the force
acting on it. This project uses the velocityVerlet algorithm
for numerical integration:
 ith
atom velocity,
 time
increment used in the simulation.
Classical potentials replace the complex quantum interaction between atoms and include experimentally fitted parameters. The reliability of MD simulations depends, mainly, on the interatomic potential chosen to describe the material. Stillinger and Weber (SW), suggested a cluster potential for modeling silicon. The pairwise potential, , (interaction between atom i, at , and atom j, at ) has the following form: (10) , A, B, p, q and a are fitting parameters. The threebody part of the SW potential (interaction between atom i, at , atom j, at , and atom k, at ) has the form: (11) where is the angle between atoms j and k subtended at vertex i. Ballamane et al. [5] suggested to modify the potential fitting parameters for silicon in order to set the minimum at 4.63eV/atom, which is the experimental value. In order to predict the brittle behavior of silicon Holland and Marder [6] proposed to stabilize the tetrahedral structure by doubling the penalty function, i.e. . The modified parameters can be found in the following file si.sw. .
Cleavage Experiment The following project performs MD computer “experiments” of dynamic crack propagation in silicon crystal using LAMMPS. The modified StillingerWeber (SW) potential is used to describe a model of a brittle crystal with a diamond structure unit cell. A microcanonical MD simulation is conducted while the sample is under constant axial strain (y axis). The initial configuration of the sample is prepared elsewhere using a code previously written by Dr. Fouad Atrash. The sample is prestrained with a 5% strain vertical to the crack propagation direction (y axis). An atomistic
striplike specimen (163840 atoms) with approximate dimensions of
1227x256x10 (in Å) is modeled. The crack is aimed to propagate on the (110) cleavage plane along
the [110] direction. A precrack, 20%
of the length of the specimen (nearly 40 atoms), was created by removing
one atomic layer at the mid of the sample's height.
The cleavage experiment (the simulation) is divided into five stages [7]: 1. sample preparation, 2. quasi static loading, 3. crack initiation modeled by MD, 4. homogenous rescaling and 5. crack propagation modeled by MD. This specimen configuration allows crack propagation
under constant driving force,
, and thus crack propagates at a constant crack velocity. The quasistatic SERR of this specimen is
calculated by the following relationship
this specimen the SERR is determined by the displacement of the boundaries, and don’t depend on the crack length. And thus the crack speed can be tuned by applying different values of .
Fig. 3 Schematics of strip like specimen Kinetic energy vs. time
The variation of
as a function of
time is shown in figure 4.
variates approximately linearly at low crack speeds and increases
wavy fluctuations around the linear trend at high speeds. This wavy
pattern is attributed to the elastic waves emitted due to crack
dynamics. The fine fluctuations lay on the wavy structure is
attributed to the thermal fluctuation of the atoms due to in initial
temperature and the thermal phonons emitted by the crack tip motion.
References: [4] F. Stillinger and T.Weber, "Computer simulation of local order in condensed phases of silicon". PRL B (1985) 52625271. [5] H. T. Balamane et al. "Comparative study of silicon empirical inetratomic potentials". PRL B 46. 22502279. [6] D. Holland and M. Marder, "Erratum (Ideal brittle fracture of silicon studied with molecular dynamics". PRL 81 (4029). [7] http://www.graduate.technion.ac.il/Theses/Abstracts.asp?Id=24738
