banner MD simulation of crack propagation in brittle crystals

118094 - Introduction to Computational Physics - Fall 2012/13
Liron Ben-Bashat 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 Born-Oppenheimer adiabatic approximation. Considering a Newtonian behavior, the equation of motion of a system consisting N particles which apply on each other a force field, Force, or equivalently potential field U is given by:

Forcefield       (8)

r, a, and m are the position, acceleration, and the mass of the i-th atom respectively and Forcei is the force acting on it.

This project uses the velocity-Verlet algorithm for numerical integration:

verlet  (9)

v - i-th atom velocity, dt - time increment used in the simulation.

Stillinger-Weber interatomic potential

Due to the absence of plastictity at the silicon crack tip, it is commonly used for atomic scale modelling of fracture.  The Stillinger-Weber interatomic classical potential [4] describes fairly well many properties of silicon but is not fitted to its elastic constants and fails to describe its brittle fracture.
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 pair-wise potential, V2, (interaction between atom i, at ri, and atom j, at rj) has the following form:

V2       (10)

sigma , A, B, p, q and a are fitting parameters.
The three-body part of the SW potential (interaction between atom i, at ri, atom j, at rj, and atom k, at rk ) has the form:

v3 (11)                             

where teta 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. lambda. 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 Stillinger-Weber (SW) potential is used to describe a model of a brittle crystal with a diamond structure unit cell. A micro-canonical 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 pre-strained with a 5% strain vertical to the crack propagation direction (y axis).

An atomistic strip-like 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 pre-crack, 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. Periodic boundary conditions are used along the crack front (z axis).

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, G0, and thus crack propagates at a constant crack velocity.

The quasi-static SERR of this specimen is calculated by the following relationship:

G0Const  (12)

eps - applied strain, H - height of the specimen and E*- effective elastic modulus. It is noted that in

 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 eps.


                                       Fig. 3 Schematics of strip like specimen

                                                 Kinetic energy vs. time

The variation of Ek as a function of time is shown in figure 4. Ek 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.


                         Fig. 4 The instantaneous average kinetic energy per atom versus time.

[4] F. Stillinger and T.Weber, "Computer simulation of local order in condensed phases of silicon". PRL B (1985) 5262-5271.
[5] H. T. Balamane et al. "Comparative study of silicon empirical inetratomic potentials". PRL B 46. 2250-2279.
[6] D. Holland and M. Marder, "Erratum (Ideal brittle fracture of silicon studied with molecular dynamics". PRL 81 (4029).