University of Southern California Dissertations and Theses

Large-scale molecular dynamics simulations of nano-structured materials

(USC Thesis Other)

Large-scale molecular dynamics simulations of nano-structured materials

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content LARGE-SCALE MOLECULAR DYNAMICS SIMULATIONS OF
NANO-STRUCTURED MATERIALS

By
PANKAJ RAJAK
_____________________________________________________________________

A Dissertation Presented to the
FACULTY OF THE USC GRADUATE SCHOOL
UNIVERSITY OF SOUTHERN CALIFORNIA
In Partial Fulfillment of the
Requirements for the Degree
DOCTOR OF PHILOSOPHY
(MATERIALS SCIENCE)

December 2018

i
Acknowledgements
This thesis dissertation brings an end to a memorable journey that has not only helped me to
become a better researcher but also helped me grew as a person. For this, I would like to express
my sincere gratitude to my advisor, Professor Priya Vashishta, for his guidance and support in this
crucial phase of my professional life. He has always stood beside me and encouraged me during
my rough times. He has entrusted me with crucial projects, for which I am immensely thankful.
He has been a wonderful advisor. I also like to express my thanks to Professor Rajiv Kalia and
Professor Aiichiro Nakano for their valuable suggestions and guidance. They persuaded me to
think more deeply about my research problems and were always available when I needed their
help.
I am also immensely grateful to Professor Ken-ichi Nomura for teaching me molecular
dynamics simulation and helping me with so many problems I had in my projects. I have cherished
the countless discussions I have had with him as they were extremely helpful. Additionally, I am
thankful to Dr. Adarsh Shekhar for his mentorship and supervision during my initial days at USC.
I am also obliged to Dr. Amit Choubey, Dr. Zaoshi Yuan, Dr. Manaschai Kunaseth and Dr. Ying
Li, who together guided me in my initial projects when I first joined the CACS group. I would also
like to thank my current colleagues at USC, Subodh Tiwari, Dr. Arvind Krishnamoorthy, Dr.
Sungwook Hong, Guoqing Zhou and Ben Watson for all for their help and support.
I owe immense gratitude to Ms. Patricia Wong, too, for taking care of all the logistics and
for being very patient with me. She was extremely helpful and understanding. The USC HPC
center and the Argonne National Lab, where I ran all my simulations, I am also thankful to their
support team as they had resolved numerous issues related to my programs. My journey at USC
would not have been a joyful ride without the help and support of my friends, thus, I would like to
ii
take this opportunity to thank my friends, Dr. Yamini Jangir, Sananda Mukherjee, Suparna
Dawalkar, Dr. Om Patri, Dr. Amulya Yadav, Swarnabha Chattaraj, Mallika Sanyal, and Bhamati
Dash for all the wonderful time we have spent together at USC.
I would also like to thank my undergrad advisor Professor Nirupam Chakraborti, at Indian
Institute of Technology, under whom I did my first research work at Kharagpur as he is the person
who sparked my interest in research. Last, but not the least, I would to express my sincerest
gratitude to my parents and my brothers, who are integral parts of my life and to whom I owe
everything.

iii
Table of Content
Acknowledgements ........................................................................................................................................................ i
List of Figures .............................................................................................................................................................. vi
Abstract ......................................................................................................................................................................... 1
Chapter 1: Introduction ............................................................................................................................................. 3
1.1 Mechanical metamaterials .......................................................................................................................... 3
1.2 Self-healing ceramics .................................................................................................................................. 4
1.3 Strain induced phase transformation in layered materials .......................................................................... 5
Chapter 2: Molecular Dynamics ................................................................................................................................ 7
2.1 Molecular dynamics overview .................................................................................................................... 7
2.2 Interaction potential .................................................................................................................................... 8
2.3 Velocity-verlet algorithm .......................................................................................................................... 11
2.4 Periodic boundary condition ..................................................................................................................... 12
2.5 Minimum image convention ..................................................................................................................... 14
2.6 Shifted potential ........................................................................................................................................ 15
2.7 Parallel molecular dynamics ..................................................................................................................... 16
2.8 Linked-list force computation ................................................................................................................... 20
Chapter 3: Mechanical Metamaterial ..................................................................................................................... 22
3.1 Background ............................................................................................................................................... 22
3.2 Simulation details...................................................................................................................................... 34
3.2.1 EAM interaction potential .................................................................................................................... 34
3.2.2 Defect identification in FCC crystal .................................................................................................... 36
3.3 Uniaxial compression of nickel nanorod and nanotube ............................................................................ 40
3.3.1 Simulation setup ................................................................................................................................... 40
3.3.2 Stress-strain curve ................................................................................................................................ 41
iv
3.3.3 Mechanical collapse of nickel nanorod ................................................................................................ 43
3.3.4 Mechanical collapse of nickel nanotube .............................................................................................. 45
3.3.5 Comparison of mechanical collapse in nanorod vs. nanotube ............................................................. 48
3.3.6 Summary .............................................................................................................................................. 51
3.4 Flat punch compression of nickel kagomé lattice ..................................................................................... 51
3.4.1 System setup......................................................................................................................................... 51
3.4.2 Stress-strain curve ................................................................................................................................ 54
3.4.3 Solid and hollow kagomé cell deformation at yield point ................................................................... 55
3.4.4 Local deformation analysis .................................................................................................................. 58
3.4.5 Solid and hollow kagomé cell deformation at higher strain ................................................................ 60
3.4.6 Comparison of kagomé cell deformation with nanorod deformation .................................................. 62
3.4.7 Bending of solid and hollow kagomé cell ............................................................................................ 63
3.4.8 Uniaxial compression of hollow kagomé lattice .................................................................................. 64
3.4.9 Summary .............................................................................................................................................. 65
3.5 Conclusion ................................................................................................................................................ 65
Chapter 4: Self-Healing Ceramics ........................................................................................................................... 67
4.1 Background ............................................................................................................................................... 67
4.2 Simulation details...................................................................................................................................... 80
4.2.1 Crystal structure of Al 2O 3 and SiC ....................................................................................................... 80
4.2.2 Interaction potential for Al 2O 3, SiC and SiO 2 ...................................................................................... 81
4.2.3 Simulation setup ................................................................................................................................... 86
4.3 Results ....................................................................................................................................................... 87
4.3.1 Stress-strain curve ................................................................................................................................ 87
4.3.2 Faceting of alumina cavity ................................................................................................................... 89
4.3.3 Crack-healing mechanism .................................................................................................................... 91
4.3.4 Secondary grain nucleation in strained Al 2O 3 ...................................................................................... 94
4.3.5 Effect of system size on crack healing ................................................................................................. 98
4.3.6 Scalable algorithm to determine fracture surface............................................................................... 100
v
4.3.7 Evaluation of mechanical property of Al 2O 3 after crack-healing by nanoindentation ....................... 100
4.4 Conclusion .............................................................................................................................................. 102
Chapter 5: Strain Induced Phase Transformation in Layered Materials ......................................................... 104
5.1 Background ............................................................................................................................................. 104
5.2 Simulation details.................................................................................................................................... 110
5.2.1 Crystal structure of MoSe 2 and WSe 2 ................................................................................................ 110
5.2.2 Stillinger-Weber interaction potential for MoWSe 2 .......................................................................... 111
5.2.3 Simulation setup ................................................................................................................................. 114
5.3 Mode-1 fracture of MoWSe 2 .................................................................................................................. 115
5.3.1 Fracture toughness of MoSe 2 and WSe 2 ............................................................................................ 115
5.3.2 Mode-1 fracture of MoWSe 2 along (10) direction ............................................................................. 117
5.3.3 Mode 1 fracture of MoWSe 2 along (11) direction ............................................................................. 118
5.3.4 Mechanism of strain induced 2H to 1T’ phase transformation .......................................................... 120
5.3.5 Mechanism of crack branching and crack closure ............................................................................. 122
5.3.6 Toughing mechanism in MoWSe 2 hetero-structure ........................................................................... 124
5.3.7 Summary ............................................................................................................................................ 125
5.4 Neural network analysis of dynamic fracture in a layered material ....................................................... 125
5.4.1 Introduction ........................................................................................................................................ 125
5.4.2 Model building ................................................................................................................................... 127
5.4.3 Visualization of the learned data representation of NN model .......................................................... 130
5.4.4 Neural network training and feature learning .................................................................................... 132
5.4.5 Hyper-parameter tuning and data visualization ................................................................................. 135
5.4.6 Summary ............................................................................................................................................ 136
5.5 Conclusion .............................................................................................................................................. 137
Conclusion ................................................................................................................................................................ 138
References ................................................................................................................................................................ 141

vi
List of Figures
Figure 2.1: Schematic of atomic trajectory simulated by MD method. ......................................................................... 7
Figure 2.2: Schematic of (a) periodic boundary condition and (b) minimum image convention. Here, red box is the
simulated MD cell, and the white boxes are the logically implemented imaginary cells. Highlighted yellow
boxes in (b) shows the regions encompassing the atoms with which the central (magenta and green) atom
interact. ............................................................................................................................................................... 13
Figure 2.3: Actual and shifted interaction potential. .................................................................................................... 15
Figure 2.4: (a) Spatial domain decomposition of PMD simulation box. (b) Shows the assigned sub-system to the
processors along the three directions. ................................................................................................................. 17
Figure 2.5: Here, each cell corresponds to a sub-system and are assigned to different processors in parallel molecular
dynamics (PMD). (a) Shows the catching layer regions (yellow) from the neighboring processors for the
highlighted red region. (b) Green arrows show migration of atoms from one processor to another processor. 18
Figure 2.6: (a) Schematic showing communication deadlock. (b) Two-phase message passing scheme to avoid
deadlock. ............................................................................................................................................................. 19
Figure 2.7: Schematic of linked list cell for force computation. .................................................................................. 20

Figure 3.1: Shows the slice of the materials-property chart for all known materials. ................................................. 23
Figure 3.2: 2D-frames showing Maxwell criteria (a) M<0, (b) M=0, (c) M>0 and (d) M=0. ..................................... 26
Figure 3.3: Dependence of the property of the lattice materials on different parameters. ........................................... 26
Figure 3.4: Comparison of modulus-density and strength-density chart in a stretching and bending dominated
structures............................................................................................................................................................. 28
Figure 3.5: Stress-strain curve for (a) a bending dominated and (b) a stretching dominated materials. ..................... 29
Figure 3.6: Fabrication of micro-lattice by self-propagating polymer waveguide.
4
.................................................... 30
Figure 3.7: Fabrication of micro-lattice by two-photon direct laser writing lithography.
34
........................................ 32
Figure 3.8: Schematic of CNA, CSP and CNP parameters for the central brown atom. Here, CNA of the central atom
in (a) is (1421), (b) opposite neighbors of the central atom are shown in same color in for CSP and (c)
common neighbor of the neighbors of the central atom is shown in yellow for CNP.
40
.................................... 38
Figure 3.9: Schematic representation of homogenous and flat punch compression. ................................................... 40
vii
Figure 3.10: Stress-Strain Curve for homogeneous and flat punch deformation of nanorod and nanotube of aspect
ratio 1:4. .............................................................................................................................................................. 42
Figure 3.11: Systematic representation of dislocation motion in nanorod. Shockley partial dislocations are shown in
yellow and hirth lock is shown in red. Various {111} planes on which dislocation motion is happening are
shown in blue, orange and green color. .............................................................................................................. 43
Figure 3.12: Systematic representation of dislocation motion nanotube. Shockley partial dislocations are shown in
red, yellow, green and magenta. ......................................................................................................................... 45
Figure 3.13: Systematic representation of dislocation motion in nanorod (a-c) and in nanotube (d-f) during
compression at mechanical collapse. .................................................................................................................. 47
Figure 3.14: Dislocation density and stress drop in the system as a function of time during the mechanical collapse
of nanorod and nanotube in homogeneous compression. ................................................................................... 49
Figure 3.15: Formation of twin region inside the deformed regions of nanotube and nanorod due to the dislocation
motion. ................................................................................................................................................................ 49
Figure 3.16: Schematic of flat punch compression of KUC. ....................................................................................... 53
Figure 3.17: Schematic of 3 × 2 × 1 kagomé lattice (KL). .................................................................................... 54
Figure 3.18: Stress-Strain of flat punch compression of SKUC and HKUC. .............................................................. 55
Figure 3.19: Plastic deformation of SKUC at yield point. ........................................................................................... 56
Figure 3.20: Dislocation density inside the (a) SKUC and (b) HKUC during yielding. ............................................. 57
Figure 3.21: Shows the regions of stress distribution with stress greater then 3GPa inside the (a) solid kagomé cell
and (b) hollow kagomé cell just before yielding. ............................................................................................... 58
Figure 3.22: Local deformation analysis of (a) SKUC and (b) HKUC showing regions of plastic deformation after
yielding. .............................................................................................................................................................. 60
Figure 3.23: Local deformation analysis showing regions of plastic deformation in SKUC at 12% and 19% strain. (a)
and (c) shows internal regions with D
2
min greater than 0.25, and (b) and (d) shows outer surface of SKUC,
which shows twin boundaries. ............................................................................................................................ 61
Figure 3.24: Local deformation analysis showing regions of plastic deformation in HKUC as a function of applied
strain. .................................................................................................................................................................. 62
Figure 3.25: Dislocation density inside the SKUC and HKUC as a function of applied strain................................... 63
viii
Figure 3.26: Angle θ between the struts at the node as a function of applied strain in SKUC and HKUC. ............... 64
Figure 3.27: (a) Stress-Strain curve of kagomé lattice, (b) regions of plastic deformation inside the KL at 15% strain.
Here, green areas are stacking fault and red lines are dislocations. ................................................................... 65

Figure 4.1: Self-healing phenomena in polymers (a) microencapsulated based healing, (b) hollow-pipes containing
healing agent and (c) microvascular network-based crack healing
55
................................................................. 71
Figure 4.2: Crack healing due to flow of liquid silica into the surface cracks. ............................................................ 76
Figure 4.3: (a) Indentation crack on Al 2O 3 surface before crack healing treatment. (b) Crack healed sample showing
fracture at the different location than the original indentation crack. ................................................................ 79
Figure 4.4: −Al 2O 3 crystal along the (a) prism plane and (b) basal plane. ................................................................ 80
Figure 4.5: (a) Schamatic of mode 1 fracture of Al 2O 3 matrix (light blue) containing a row of n-SiC (green) with
SiO 2 shells (red) with a pre-crack (dark blue). (b) Shows the boundary conditions in different regions of the
system. ................................................................................................................................................................ 87
Figure 4.6: (a) Effect of inter-particle distance on crack healing Al 2O 3 composite with single row of of SiC/SiO2
nanoparticles. (b) Stress-Strain curve for Al 2O 3 composite with four rows of SiC/SiO2 nanoparticles and for
the pure Al 2O 3 system during mode 1 fracture. Pure Al 2O 3 and systems with larger inter-particle distance
fractures at 3.5% however, alumina composite with nanoparticles withstand up to 5% strain.......................... 88
Figure 4.7: Faceting of an Al 2O 3 cavity around SiC/SiO 2 nanoparticle at 3.4% strain. Here, Al 2O 3 contains two rows
of SiC/SiO 2 nanoparticle. (a) Al 2O 3 cavity further away from the pre-crack changes to hexagonal shape while
(b) Al 2O 3 cavity near the pre-crack changes from spherical to pentagonal shape. These facets form along
prismatic (A) 2110 and prismatic (M) 1010 planes. ......................................................................................... 90
Figure 4.8: (a) Crack healing in Al 2O 3 by flow of molten silica. Here, the alumina matrix is black and yellow, and
red and green spheres are Si, O and C atoms, respectively. The NP diameters is 7 nm and Al 2O 3 contains 1,2,3
and 4 rows of NPs. Each row contains 2 NPs of diameter 7 nm seperated by 2 nm. (b) Self-diffusion
coefficients of Si and O as a function of the applied strain. (c) Strain dependence of Δσ xx, which is the
difference in the stress at the crack tip and NPs. ................................................................................................ 91
ix
Figure 4.9: (a) Distance between crack-front and NPs and (b) distance between NPs dp-p in an Al 2O 3 matrix
containing a row of two nSiC/SiO 2 nanoparticle. The initial separation between the crack and nanoparticles is
4.5 nm, and the diameter of each nanoparticle is 7 nm. ..................................................................................... 93
Figure 4.10: Nucleation and growth of multiple secondary grains at the interface of an Al 2O 3 –nSiC/SiO 2. Here, the
Al 2O 3 matrix contains two rows of nSiC/SiO 2 nanoparticles. Secondary grains form above 3% strain and
facets of the Al 2O 3 cavity act as nucleation sites for these grains. Particle 1 in the figure is in row 1, which is
closer to the crack then particle 2 in row 2. ........................................................................................................ 94
Figure 4.11: Analysis of stress distribution and local damage inside the Al 2O 3 matrix containing 2 rows of NPs at
3.4% strain. (a) Stress distribution is plotted on a plane at z=20 nm. (b) Local damage on a cross-section of
Al 2O 3 at z=20 nm. Here, yellow and red are silicon and oxygen atoms, respectively. Al 2O 3 atoms with D
2
min
value less than 0.25 are shown in black and the remaining Al 2O 3 atoms are colored by their D
2
min value. ...... 95
Figure 4.12: (a) Local deformation analysis of the deformed region in Al 2O 3 matrix at 4.5% strain. (b) Comparison
of Al-O-Al and O-Al-O bond angle distributions in the deformed region of alumina. (c) Comparison of radial
distribution function g(r) for Al-O and coordination number n(r) in the deformed region. Here, dotted lines are
g(r) and n(r) plots in Al 2O 3 matrix at 0% strain and solid lines in the deformed region of Al 2O 3 matrix at 4.5%
strain. .................................................................................................................................................................. 98
Figure 4.13: Effect of (a) 2 rows and (b) 4 rows of NPs on crack healing in Al 2O 3. Here, Al 2O 3 atoms are shown in
black color, and Si and O atoms of NPs are shown in yellow and red, respectively. Atoms near the crack
surface or damaged zone are colored by their D
2
min value. Crack formation happens due to void nucleation and
coalescence along the slip direction [0111] of the prism plane of -Al 2O 3, as shown in Figure 4.12(a). ........ 99
Figure 4.14: Load-displacement curve for the nanoindentation of Al 2O 3 before and after crack healing at room
temperature. ...................................................................................................................................................... 102

Figure 5.1: Synthesis of (a) MoS 2 and (b) MoWSe 2 monolayer by CVD.
86,89
.......................................................... 106
Figure 5.2: (a) Schematic of the mechanical straining of MoWSe 2 monolayer by three-point bending experiment. (b)
Shift in Raman mode of MoSe 2 and WSe 2 as a function of applied strain....................................................... 107
Figure 5.3: Schematic of the atomic displacement of four active Raman-modes in MoS 2.
90
.................................... 108
Figure 5.4: (a) Schematic of mode 1 fracture in MoWSe 2 alloy containing a single patch of WSe 2. (b) and (c) shows
atomic view when loading direction is perpendicular to (1,1) and (1,0), respectively. ................................... 115
x
Figure 5.5: K IC value along different loading direction in MoSe 2 and WSe 2 ............................................................ 116
Figure 5.6: Distribution of the stress component  yy during crack propagation in MoWSe 2 monolayer for (a) (1,0)
and (b) (1,1) loading direction. Critical stresses near the crack tips reach 20 GPa at the onset of fracture. .... 117
Figure 5.7: Stress induced 2H-1T phase transformation near the tip of the propagating crack. (a) Snapshot of the 2H-
1T phase transformation as the crack propagates through the MoSe 2/WSe 2 alloy. The highlighted region
indicates the two different crystal structures in the process zone. Green atoms are at the interface of the 2H
phase and the disordered 1T phase. (b) Simulated Raman spectra of the MoSe 2/WSe 2 alloy during crack
propagation. Spectra obtained during the initial stages of crack propagation (0-1 ps) are characteristic of the
2H crystal structure of the MoSe 2/WSe 2 alloy, while spectra from later times are obtained from simulation
cells with greater 1T content. Figures (c)-(e) are snapshots from the MD simulation of atomic configurations
in the process zone. These snapshots show atomic displacements corresponding to the motion of the 2H-1T
phase transformation front (dashed brown line). Here Mo, Se, and W are shown in pink, yellow and blue
colors, respectively. The transformation of Mo atoms (1, 4, 7) to the 1T structure is due to the displacement of
Se atoms (2, 3, 5, 6, 8 and 9) highlighted in figures (c)-(e). Dotted blue curve indicates bond breaking and
dotted blue lines shows bond formation. For example, in figure (c) the formation of 1(Mo)-5(Se) bond and
breaking of 1(Mo)-2(Se) and 7(Mo)-5(Se) bonds transforms 1(Mo, cyan) from 2H to 1T structure. The
resultant 1T structure is shown in figure (d). During this transformation, the Se pairs 2-3 and 5-6 move in
opposite direction as shown in figure (d). ........................................................................................................ 119
Figure 5.8: Local deformation around the crack tip in a strained in an MoSe 2/WSe 2 monolayer in [1,1] configuration.
(a) Local strains in the process zone around the pre-crack and crack branches. In an unstrained MoWSe 2
monolayer in the 2H or 1T phase, Mo and W atoms belong to six and Se to three rings, each with four
members. Here, green dots represent defects, i.e., Mo and Se atoms in which one of the 4-membered ring is
replaced by a larger ring. .................................................................................................................................. 121
Figure 5.9: (a-d) shows crack propagation and crack branching mechanism in MoSe 2/WSe 2 monolayer. (a) At the
crack tip, multiple W (orange)-Se (cyan), W (orange)-Se (light blue) and W (pink)-Se (light blue) bonds
breaking leads to the formation of defects at the crack tip along different directions as shown in (b) which
eventually leads into (c-d) crack branching in the material.............................................................................. 123
xi
Figure 5.10: Area of the process zone in the MD simulation as a function of time in a pure MoSe 2 monolayer (blue)
and in a MoWSe 2 alloy (red). ........................................................................................................................... 124
Figure 5.11: Schematic of the neural network model showing the mapping of atomic coordinates into sets of radial
and angular feature vectors, NN layer architecture and the six out phases predicted by the model. ............... 129
Figure 5.12: (a-c) shows PCA and (d-e) shows t-SNE analysis of the input feature vector and the learned
representation of data in hidden layer 1 and 3 of the neural network. ............................................................. 130
Figure 5.13: (a) Shows 10 features whose gradient change is maximum during training using full feature vector set.
(b) 10 features whose gradient change is minimum during training. (c) and (d) shows training and validation
accuracy of models, which are training with 50 features with maximum gradient change or minimum gradient
change in the original feature space. ................................................................................................................ 133
Figure 5.14: 10 most significant radial and angular symmetry functions learned by the neural network. ................ 135
Figure 5.15: Visualization of validation set of a MoWSe 2 fracture data for modes tuned with two different cutoffs for
angular symmetry function calculation (a) r cut= 8 Å and (b) r cut=3 Å. .............................................................. 136

1
Abstract
Nanostructured materials have exceptional mechanical properties and have a wide variety of
applications ranging from manufacturing, aerospace to healthcare. However, the service life of
these materials is often limited due to the formation of cracks and defects that leads to their
eventual failure, thus insights into atomistic mechanisms underlying these processes are valuable
inputs for the design of next-generation nanostructured materials. In this research work, large-scale
molecular dynamics simulation, combined with machine learning and data mining techniques is
used to understand three basic problems related to nanostructured materials design with superior
mechanical properties, which are: (1) ultra-light material design using hierarchical structures like
kagomé lattices, (2) self-healing of cracks in alumina composite embedded with SiC/SiO2
nanoparticles, and (3) toughing of MoWSe2 monolayer during crack propagation by novel strain-
induced structural phase transformation.
Many natural and man-made materials exhibit hierarchical structures like for example bird
wings, Eiffel tower, Garabit Viaduct. This hierarchy in structures can be used for materials design
at the atomic level as it allows us to tune several materials properties like stiffness, compressive
strength, fracture toughness by just changing the underlying structural arrangement of materials.
In the first work, we have studied the deformation behavior of an ultralight architecture consisting
of hollow Ni nanotubes and solid nanorods arranged as a 3-D kagomé structure. Our result shows
that hollow kagomé structure have higher strength and stiffness in comparison of solid kagomé
structures. The hollow kagomé structure shows two stage deformation mechanism where up to
11% strain, deformation is highly localized near the nodal region and at higher strain they show
bending-dominated deformation. Whereas in solid kagomé structures, large defamations are
observed in the entire structure even at lower strain.
2
Designing material systems that can sense and repair damage is of great interest. Inspired by
biological systems, material scientists are developing self-healing materials in which damage
initiation triggers an automatic release and transport of healing agent to the damage site and repairs
it. In the second work, we have investigated crack healing in a ceramic nanocomposite of alumina
containing SiC/SiO2 nanoparticles. Our result shows that the nanoparticles closest to the advancing
crack in alumina distort to create nanochannels through which silica flows toward the crack and
heals the damage. Under applied strain, alumina matrix also forms facets at interfaces. These facets
act as nucleation sites for multiple secondary amorphous grains in alumina matrix. Nanovoids form
along the grain boundaries but these nanovoids are again healed by molten silica released from the
nanoparticles. The mechanisms of crack healing in alumina will help experimentalists to design
novel self-healing ceramics for a broad class of energy technologies
Two-dimensional materials exhibit different mechanical property in comparison of their bulk
counterpart owing to their monolayer atomic thickness. In fact, properties of these layered
materials can be tuned by mechanical deformation. In the third work, we have examined the
atomistic mechanisms of defect formation and crack propagation in a strained monolayer MoWSe2
alloy using molecular dynamics. Our result shows that under strain, the crack meanders as it
propagates through the system and branches into daughter cracks when it crosses the MoSe2/WSe2
interface and enters the WSe2 patches. The most dramatic change occurs in the process zone around
the crack tip, where large stress concentration and the resulting biaxial tensile strain trigger an
irreversible local structural phase transformation from the ground state 2H crystal structure to the
1T crystal structure. This phenomenon of crack branching, crack closure and strain induced phase
transformation in MoWSe2 alloy increases its fracture toughness.

3
Chapter 1: Introduction
1.1 Mechanical metamaterials
Hierarchy in structures exists all around us. Many natural and man-made materials exhibit
this hierarchy. A common example of hierarch
1
ical structure found in nature is birds’ nests.
2,3
At
the macro level, bird nest is formed by leaves of grass, which are wavy and bend when nest
deforms. And at the finer level, these grass leaves are composed of three-dimensional form like
arrangement of cells with wavy sides. As a result of which the cell walls of the grass leaf bends
too under strain and make bird nest very damage tolerant under storms.
Taking cue of hierarchical structural design of nature, material scientists are designing
micro-architecture material with hierarchical structures. These micro-architecture materials are
able to occupy regions of the materials property space which are previously been un-explored.
4-8

We are now able to design materials with properties like high-strength, high stiffness and high
fracture toughness at low density. These micro-architectured materials, which are based on the
scale of hierarchy and nodal connectivity can be tuned to provide wide range of materials with
desired stiffness, strength and fracture toughness.
2
This hierarchy of structures can also be used to
design materials with negative poisson’s ratios and super-plasticity.
Mechanical properties of these micro-architectured materials are dependent upon many
factors which are the constituent base material, structure of the micro-lattice, density, wall
thickness, etc. For example, a cellular lattice of ceramic material made from octet truss is highly
ductile when the wall thickness of struts is below certain limit. The effects of these design
parameters on the underlying atomistic deformation mechanism of micro-architectured materials
is not well understood. Hence, here we have used molecular dynamics simulation to understand
the effect of wall thickness on the atomistic deformation mechanism of nickel kagome lattices.
4
1.2 Self-healing ceramics
There is a great deal of interest in designing material systems having the ability to sense and
repair damage.
1,9
Usually damage initiation involves the formation of defects such as dislocations,
voids and micro-cracks. These defects grow and coalesce to cause material failure. If we can detect
and repair these micro-cracks at early stage, then it will enhance their service life, which in many
situations is not a feasible solution because it requires constant human intervention. In this respect,
we can learn a lot from nature as many natural materials shows automatic healing phenomena, for
example, plants and animals have the ability to heal their wound. Biological materials have two
unique characteristics that helps them to achieve high strength and damage tolerance. These are
hierarchy in their micro-structure and self-healing phenomena.
Taking this cue from biology, material scientists are developing novel self-healing systems
to enhance reliability and lifetime of materials while reducing the cost of manufacturing,
monitoring and maintenance.
10,11
Most of the approaches to designing self-healing material
systems involve nanoscale or mesoscale containers filled with a healing material, which is released
upon damage initiation either by an external stimulus (heat, light, electric or magnetic field) or an
internal change in the state of the system (pH, stress or temperature).
11-14
The healing material is
transported to the damage site and immobilized after healing damage. Different materials require
different self-healing conditions — ambient conditions for concrete, low temperatures for
polymers, relatively high temperatures for metals, and very high temperatures for ceramics.
In comparison of polymers, self-healing phenomena in ceramics is harder and is not
extensively studied because they operate at high temperature, have high melting point and hence
low atomic mobility.
15
Crack healing in ceramics is usually done by oxidation. All ceramics that
shows oxidation-based crack healing, contains two ceramic phases where one phase reacts with
5
the oxygen from the air and forms oxide. Examples of crack healing by oxidation in ceramics
includes mullite-SiC composite, Si3N4-SiC composite,
16,17
Al2O3-SiC composite
18,19
and MAX
phase ceramics (Ti3AlC2, Ti2AlC).
20,21
However, very little is known about the mechanism of crack
healing at atomistic level for example the effect of nano-particle size and distribution inside the
ceramic matrix. Here, we have studied the crack healing mechanism of Al2O3 composite containing
SiO2 coated SiC nanoparticles using molecular dynamics simulation.
1.3 Strain induced phase transformation in layered materials
Since the advent of single layer graphene, a wide variety of two-dimensional (2D)
materials, especially semiconducting transition metal dichalcogenides (TMDCs), have been
isolated with a suite of interesting properties and applications.
22-27
These TMDCs have a range of
band gap, for example, the semiconductors with large band gap (MoS2 ~ 2 eV) and
superconducting PDTe2. TMDCs have application in many areas which are catalyst in
dehydrosulfurization, dry lubricants, photovoltaic power generation and as an ingredient in some
Li-ion batteries.
TMDC monolayers have interesting mechanical properties as well. For example, Young’s
modulus of MoS2 and WS2 monolayer calculated from the force-displacement curve of
nanoindentation experiment is shown to be ~270 GPa, which is higher than the Young’s modulus
of bulk MoS2 (~240 GPa) and WS2 (~170 GPa). Several efforts have been made recently to tune
the optical bandgap of the atomically thin sheets of TMDCs by mechanical deformation. Similarly,
the alloyed system of the TMCD monolayers have interesting mechanical properties. Examples of
such materials include the simultaneous presence of Mo and W chalcogenides in a monolayer.
However, the strain effects on the monolayer TMDC alloys have not been explored well due to
6
the experimental challenges resulted from the high sensitivity of TMDCs under electron beam and
limited deformation in atomic force microscopy (AFM).
Here, we have used molecular dynamics simulation to understand the mechanism of crack
propagation in MoWSe2 monolayer for different loading directions. To analyze the fracture
simulation results, a neural network model is built, which automatically learns the relevant
information of the simulation data like phase transformation, defect formation.

7
Chapter 2: Molecular Dynamics
2.1 Molecular dynamics overview
Molecular dynamics (MD) is a computer simulation technique which is used to study the
physical and chemical processes of materials at the atomic resolution. In MD, we simulate the
trajectory of atoms/molecules as a function of time by solving Newton’s equation of motion that
are second order non-linear coupled ordinary differential equations. Using the trajectory
information, we can study various physical phenomena, for example, brittle vs. ductile fracture,
tensile vs. compressive deformation, response to shock, oxidation and chemical reactions.

Figure 2.1: Schematic of atomic trajectory simulated by MD method.
As stated above, in MD we solve Newton’s equation of motion. For a system of N particles,
the positions and velocities of the N atoms, { (𝑟 𝑁 (0)
⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗
,𝑣 𝑁 (0)
⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗
)|𝑖 =0,…,𝑁 −1} are required as
initial inputs. Using this information, it computes the forces on each atom from all neighboring
atoms within prescribed cutoff distance and calculate the acceleration of atoms. Then the atomic
position and velocities are updated based on the obtained acceleration with positions and velocities
at previous time step. This completes a single iteration of MD and is repeated to simulate the
trajectory of atoms as a function of time, { (𝑟 𝑁 (0)
⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗
,𝑣 𝑁 (0)
⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗
)→(𝑟 𝑖 (𝑡 1
)
⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗
,𝑣 𝑖 (𝑡 2
)
⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗
)→⋯→
8
(𝑟 𝑖 (𝑡 )
⃗⃗⃗⃗⃗⃗⃗⃗
,𝑣 𝑖 (𝑡 )
⃗⃗⃗⃗⃗⃗⃗⃗⃗
)|𝑖 =0,…,𝑁 −1}. Table 1 shows the pseudocode of the MD simulation, where Fi(t)
and ai(t) are the force and acceleration of the i
th
atom, respectively, and m is the mass of atom.

Table 1: Pseudocode of a molecular dynamics simulation.
In Table 1, the key ingredient is the interaction potential function between atoms, which
should model underlying physics of the target system well. The obtained forces based on the
interaction potential are used to update the position and velocities of atoms with time integration
algorithms like leaf-frog, velocity-verlet, etc. A brief description of the interaction potential and
velocity-verlet algorithm is given the subsequent sections.
2.2 Interaction potential
In MD, forces between atoms are calculated by taking the derivate of the potential energy,
𝑉 (𝑟 𝑁 ⃗⃗⃗⃗ ) , with respect to atomic positions. Here, 𝑉 (𝑟 𝑁 ⃗⃗⃗⃗ ) is the function of atomic positions of N atoms.
Hence, the force the k
th
atom will be the rate of change in 𝑉 (𝑟 𝑁 ⃗⃗⃗⃗ ) when we slightly change the
coordinate of k
th
atom while keeping the coordinate of all other atoms fixed.
𝐹 𝑘 ⃗⃗⃗⃗
= −
𝜕 𝜕 𝑟 𝑘 ⃗⃗⃗⃗
𝑉 (𝑟 𝑁 ⃗⃗⃗⃗ )=−(
𝜕𝑉 (𝑟 𝑁 ⃗⃗⃗⃗ )
𝜕 𝑥 𝑘 ,
𝜕𝑉 (𝑟 𝑁 ⃗⃗⃗⃗ )
𝜕 𝑦 𝑘 ,
𝜕𝑉 (𝑟 𝑁 ⃗⃗⃗⃗ )
𝜕 𝑧 𝑘 ) Eq. 1
Where the potential energy 𝑉 (𝑟 𝑁 ⃗⃗⃗⃗ ) is
𝑉 (𝑟 𝑁 ⃗⃗⃗⃗ )=∑ 𝑢 (𝑟 𝑖𝑗
)
𝑖 <𝑗 =∑ ∑ 𝑢 (|𝑟 𝑖𝑗
⃗⃗⃗ |)
𝑁 −1
𝑗 =𝑖 +1
𝑁 −2
𝑖 =0
Eq. 2
Input: {(𝑟 𝑖 (0)
⃗⃗⃗⃗⃗⃗⃗⃗⃗
,𝑣 𝑖 (0)
⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗
)| 𝑖 =0,…,𝑁 −1
Step t = 1, N_step_max
1. Compute 𝐹 𝑖 (𝑡 )
⃗⃗⃗⃗⃗⃗⃗⃗⃗
for all atoms
2. Compute 𝑎 𝑖 (𝑡 )
⃗⃗⃗⃗⃗⃗⃗⃗⃗
for all atoms, 𝑚 𝑎 𝑖 (𝑡 )=𝑚 𝑟 𝑖 (𝑡 )
⃗⃗⃗⃗⃗⃗⃗⃗
̈ =𝐹 𝑖 (𝑡 )
⃗⃗⃗⃗⃗⃗⃗⃗⃗

3. Update position and velocity, (𝑟 𝑖 (𝑡 )
⃗⃗⃗⃗⃗⃗⃗⃗
,𝑣 𝑖 (𝑡 )
⃗⃗⃗⃗⃗⃗⃗⃗⃗
) for all atoms using 𝑎 𝑖 (𝑡 )
⃗⃗⃗⃗⃗⃗⃗⃗⃗

and (𝑟 𝑖 (𝑡 −1)
⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗
,𝑣 𝑖 (𝑡 −1)
⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗
)

9
Equation 2 shows the total potential energy of the system, 𝑉 (𝑟 𝑁 ⃗⃗⃗⃗ ) , in terms of individual
potential energy of the atoms, 𝑢 (𝑟 𝑖𝑗
) . Here, the potential energy of an atom can be described by
different types of interaction potential. For example, it can be represented as a simple pair-wise
function such as Lennard-Jones (L-J) potential, which is suitable to model materials like liquid
argon and inert gas. Other types of interaction potentials that exist are non-reactive and reactive
potentials. Non-reactive potentials can model bond stretching, bond bending and torsion but do
not allow bond breaking and formation, whereas reactive potentials as ReaxFF are designed to
describe chemical reactions well. Simple interaction potentials (L-J, Morse potential) have few
parameters (2-4) that are fitted using experimental data to get correct lattice constant, cohesive
energy and bulk modules. More advanced force fields (ReaxFF, EAM) have many tunable
parameters to accurately describe elastic constants, bond dissociation energy and atomic charges.
Depending on the physical processes we want to simulate, we choose different types of interaction
potential. For example, ReaxFF allows bond breaking and bond formation, and thus it is a desired
interatomic potential to simulate phenomena involving chemical reactions whereas EAM
potentials describes elastic properties well and are good for mechanical deformation simulation.
A brief description of two interaction potentials, L-J potential and Smith and Bharadwaj (SB)
potential (non-reactive force field), are given below.
𝑢 (𝑟 )=4𝜖 ((
𝜎 𝑟 )
12
−(
𝜎 𝑟 )
6
) Eq. 3
L-J potential (Eq. 3) has two tunable parameters, 𝜎 and 𝜖 , which are fitted to get correct
cohesive energy and lattice constant. It consists of two terms, the first term allows short range
repulsion, and the second term is for long range interaction. The Pauli exclusion principal tells that
two electrons cannot occupy the same state as it causes repulsive interaction between atoms when
10
atoms approach very close to each other. This property of atoms is described by the first term. The
second long-range attractive term describes van der Waals interaction between atoms.
Smith and Bharadwaj (SB) potential (Eq. 4) is a non-reactive force field which has been
used energetic molecular crystals such as RDX and HMX. It consists of bonding, non-bonding and
Coulombic interaction terms. Bonding interaction includes bond stretching and bond bending.
Non-bonding interaction is descried by L-J and Buckingham potentials. The Coulombic interaction
is represented by fixed atomic charges, and the long-range electrostatic interactions are calculated
using the Particle-Particle Particle Mesh (PPPM) method. The crystal density, coefficient of linear
expansion and elastic constants predicted by the SB potential agree well with experimental results.
This interatomic potential has been used extensively to study shock response and mechanical
properties of RDX and HMX. The functional form of the SB potential is shown in equation 4.
𝑈 =∑
1
2
𝐾 𝑖𝑗
𝑆 (𝑟 𝑖𝑗
−𝑟 𝑖𝑗
0
)
𝑏𝑜𝑛𝑑𝑠 2
+∑
1
2
𝐾 𝑖𝑗𝑘 𝐵 (𝜃 𝑖𝑗𝑘 −𝜃 𝑖𝑗𝑘 0
)
𝑎𝑛𝑔𝑙𝑒 2
+
∑
1
2
𝐾 𝑖𝑗𝑘𝑙 𝑇 𝑝𝑟𝑜𝑝𝑒𝑟 𝑑𝑖 ℎ𝑒𝑑𝑟𝑎𝑙𝑠 [1−𝑐𝑜𝑠 (𝜂 𝜙 𝑖𝑗𝑘𝑙 )]+∑
1
2
𝐾 𝑖𝑗𝑘𝑙 𝐷 𝛿 𝑖𝑗𝑘𝑙 2
𝑖𝑚𝑝𝑟𝑜𝑝𝑒𝑟 𝑑𝑖 ℎ𝑒𝑑𝑟𝑎𝑙𝑠 +
∑ ∑ ((𝐴 𝑖𝑗
𝑒𝑥𝑝 (−𝐵 𝑖𝑗
𝑟 𝑖𝑗
)−
𝐶 𝑖𝑗 𝑟 𝑖𝑗 6
+𝑘 𝑒 𝑞 𝑖 𝑞 𝑗 𝑟 𝑖𝑗 )
𝑁 𝑗 >𝑖 𝑁 −1
𝑖 =1
Eq. 4
The first four terms represent bonding interaction, while the last term is for non-bonding
(Buckingham potential) and Coulombic interactions. Here, rij is the distance between the i
th
and j
th

atoms, 
ijk
is the bond angle formed by the i
th
, j
th
and k
th
atoms, 
ijkl
is the dihedral angle formed
by the i
th
, j
th
, k
th
and l
th
atoms, and qi is the charge of the i
th
atom. In equation 4, 𝐾 𝑖𝑗
𝑆 , 𝐾 𝑖𝑗𝑘
𝐵 , 𝐾 𝑖𝑗𝑘𝑙 𝑇 ,
𝐾 𝑖𝑗𝑘𝑙 𝐷 are the force constants for bond-stretching, bond-bending, torsion and out-of-plane bending,
respectively; 𝑟 𝑖𝑗
0
,𝜃 𝑖𝑗𝑘
0
are equilibrium bond length and bond angle, respectively; 𝐴 𝑖𝑗
,𝐵 𝑖𝑗
,𝐶 𝑖𝑗
are
parameters of the Buckingham potential; and  and ke are constants. Since, the Buckingham
11
potential becomes attractive at small interatomic distances, the L-J interaction is further added in
the interaction potential to describe the non-bonding interaction.
2.3 Velocity-verlet algorithm
In order to compute the trajectories of atoms on a computer, we need to discretize the
atomic trajectory. Hence, instead of considering atomic trajectory as continues function in time
{(𝑟 𝑖 (𝑡 )
⃗⃗⃗⃗⃗⃗⃗⃗
,𝑣 𝑖 (𝑡 )
⃗⃗⃗⃗⃗⃗⃗⃗⃗
),𝑡 ≥0}, we consider it as a sequence of discrete states, {(𝑟 𝑖 (0)
⃗⃗⃗⃗⃗⃗⃗⃗⃗
,𝑣 𝑖 (0)
⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗
)→
(𝑟 𝑖 (Δ)
⃗⃗⃗⃗⃗⃗⃗⃗⃗
,𝑣 𝑖 (Δ)
⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗
)→(𝑟 𝑖 (2Δ)
⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗
,𝑣 𝑖 (2Δ)
⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗
)→⋯}, where the difference between the two consecutive states
is equal to the time step, Δ, of the MD simulation. One such algorithm to calculate this discrete
sequence of atomic state during MD simulation is Velocity-Verlet Algorithm. The pseudo-code of
the algorithm is given in Table 2.

Table 2: Pseudocode of a velocity-verlet algorithm.
The Velocity-Verlet algorithm takes as input the initial position and velocities of all the
atoms in the system and compute their initial acceleration ai using the function computeAccel.
ComputeAccel calculates the force on atoms from their neighbors by taking the derivate of the
potential energy, which is then used to compute ai. Once we have the initial accelerations, step 1
to 4 is repeated in a loop for a total of N
step_max
steps where in each iteration the atomic positions,
Velocity-Verlet algorithm
Input: {(𝑟 𝑖 (0)
⃗⃗⃗⃗⃗⃗⃗⃗⃗
,𝑣 𝑖 (0)
⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗
)| 𝑖 =0,…,𝑁 −1
Compute 𝑎
𝑖 as a function of {𝑟
𝑖 } for all atoms 𝑖 =0,…,𝑁 −1 using computeAccel()
Step t = 1, N_step_max
1. 𝑣
𝑖 ⟵ 𝑣
𝑖 +𝑎
𝑖 Δ
2
𝑓𝑜𝑟 𝑎𝑙𝑙 𝑎𝑡𝑜𝑚𝑠
2. 𝑟
𝑖 ⟵𝑟
𝑖 +𝑣
𝑖 Δ 𝑓𝑜𝑟 𝑎𝑙𝑙 𝑎𝑡𝑜𝑚𝑠
3. Compute 𝑎
𝑖 as a function of {𝑟
𝑖 } for all atoms using computeAccel ()
4. 𝑣
𝑖 ⟵ 𝑣
𝑖 +𝑎
𝑖 Δ
2
𝑓𝑜𝑟 𝑎𝑙𝑙 𝑎𝑡𝑜𝑚𝑠

12
velocities and accelerations are updated. Here each iteration is equivalent to Δ timestep, and hence
total simulation time will be Δ× N
step_max
. We can observe from Table 2 that the velocity is
updated twice in each iteration.
It is very important to choose the timestep, Δ, carefully. If the timestep is too small, it will
take a very long time to simulate any physical or chemical phenomenon. On the other hand, if the
time-step is too large then the system will blow up due to numerical instability. In general, the time
step chosen in MD simulations are in the order of femto-second, (~10
−15
𝑠𝑒𝑐 =1𝑓𝑠 ) . The
reasoning behind is that the time step of MD simulation should be well below the period of the
highest vibrational frequency of the atoms in the system.
2.4 Periodic boundary condition
The number of atoms that we simulate in MD is finite. Hence, to simulate bulk properties,
i.e., systems with infinite number of atoms, we need to apply some kind of boundary conditions to
our MD simulation box. Different types of boundary conditions that exists are fixed boundary
condition (FBC), shrink-wrapped boundary condition (SWBC) and periodic boundary condition
(PBC). In the fixed boundary condition, when an atom moves outside the MD simulation box then
it is lost, thus the total number of atoms may not remain constant in these kinds of simulation.
However, in shrink-wrapped boundary condition, when an atom moves outside the MD simulation
box then the shape of MD simulation box is changed so as to encompass that atom again within
the simulation cell. In both of these boundary conditions, atoms across the boundary do not interact
with each other. Hence, if an atom crosses the boundary at one side of the simulation cell then it
does not enter the simulation cell from the other side. These types of boundary conditions are
called non-periodic boundary conditions and are used to simulate surfaces and properties of
clusters.
13

Figure 2.2: Schematic of (a) periodic boundary condition and (b) minimum image convention.
Here, red box is the simulated MD cell, and the white boxes are the logically implemented
imaginary cells. Highlighted yellow boxes in (b) shows the regions encompassing the atoms with
which the central (magenta and green) atom interact.

Bulk properties are simulated using periodic boundary condition (PBC). In PBC, when an
atom crosses the MD simulation box, it enters back the simulation box from the other side. Also,
here atoms interact with each other across the boundary. Hence, PCB can simulate infinitely large
system without any surface effect on materials property. PCB represents the bulk system with
infinite number of cells where the behavior of atoms in each cell is identical. We simulate only
one these cell and rest are implemented logically through PBC condition. Figure 2.2(a) shows a
schematic of PBC where the highlighted red cell shows the area that contains real atoms and the
white cells its images. Equation 5 shows a pseudo-code implementing the PBC condition, which
states that if the coordinate of an atom becomes larger then MD simulation box size then subtract
the box size from its coordinate. Similarly, if the coordinate of an atom becomes less than zero,
we add MD simulation box length to its position. In other words, it means that a periodic image of
the atom of the simulated cell has entered the simulation box from the other side (highlighted with
blue arrow in Figure 2.2a, where the magenta atom which is an image of the yellow atom enters
the simulation from the other side when the yellow atom moves out of the box).
14
{

𝑥 𝑖 ←𝑥 𝑖 −𝑆𝑖𝑔𝑛𝑅 (
𝐿 𝑥 2
,𝑥 𝑖 )−𝑆𝑖𝑔𝑛𝑅 (
𝐿 𝑥 2
,𝑥 𝑖 −𝐿 𝑥 )
𝑦 𝑖 ←𝑦 𝑖 −𝑆𝑖𝑔𝑛𝑅 (
𝐿 𝑦 2
,𝑦 𝑖 )−𝑆𝑖𝑔𝑛𝑅 (
𝐿 𝑦 2
,𝑦 𝑖 −𝐿 𝑦 )
𝑧 𝑖 ←𝑧 𝑖 −𝑆𝑖𝑔𝑛𝑅 (
𝐿 𝑧 2
,𝑧 𝑖 )−𝑆𝑖𝑔𝑛𝑅 (
𝐿 𝑧 2
,𝑦 𝑖 −𝐿 𝑧 )
Eq. 5
Where
𝑆𝑖𝑔𝑛𝑅 (
𝐿 𝑥 2
,𝑥 𝑖 )={

𝐿 𝑥 2
𝑥 𝑖 >0
−
𝐿 𝑥 2
𝑥 𝑖 ≤0

2.5 Minimum image convention
An atom interacts with all the atoms in the simulation cell and also with their images in the
logically implemented cells. In fact, it can interact with its own images except for itself. However,
the forces on an atom from its neighbor that are at very large distance are negligible, and thus can
be ignored. Hence in practice, the size of the MD simulation box (𝐿 𝑥 ×𝐿 𝑦 ×𝐿 𝑧 ) is made large
enough so that we can ignore contributions from atoms with 𝑢 (𝑟 >
𝐿 =max(𝐿 𝑥 ×𝐿 𝑦 ×𝐿 𝑧 )
2
). It implies
that atoms in the central MD box interact at max with the nearest image cells. In other words, for
each atom in the MD simulation box, we can imagine that it is surrounded by a cube of size
𝑙 𝑥 ×𝑙 𝑦 ×𝑙 𝑧 with it at the center of the cube, and it interacts with all the atoms enclosed by this
cube (see Figure 2.2b). For atoms at the boundary of the MD simulation box, it means that it will
interact with atoms in the simulation MD box and with the atoms from the first image cell
surrounding the original simulation box (see Figure 2.2b). Mathematically, minimum image
convention is implemented as:
15
{

𝑥 𝑖𝑗
←𝑥 𝑖𝑗
−𝑆𝑖𝑔𝑛𝑅 (
𝐿 𝑥 2
,𝑥 𝑖𝑗
+
𝐿 𝑥 2
)−𝑆𝑖𝑔𝑛𝑅 (
𝐿 𝑥 2
,𝑥 𝑖𝑗
−
𝐿 𝑥 2
)
𝑦 𝑖𝑗
←𝑦 𝑖𝑗
−𝑆𝑖𝑔𝑛𝑅 (
𝐿 𝑦 2
,𝑦 𝑖𝑗
+
𝐿 𝑦 2
)−𝑆𝑖𝑔𝑛𝑅 (
𝐿 𝑦 2
,𝑦 𝑖𝑗
−
𝐿 𝑦 2
)
𝑧 𝑖𝑗
←𝑧 𝑖𝑗
−𝑆𝑖𝑔𝑛𝑅 (
𝐿 𝑧 2
,𝑧 𝑖𝑗
+
𝐿 𝑧 2
)−𝑆 𝑖𝑔𝑛𝑅 (
𝐿 𝑧 2
,𝑦 𝑖 −
𝐿 𝑧 2
)
Eq. 6
Figure 2.2b shows a schematic of the force computation on an atom at center of the MD
box (magenta atom) and at the boundary (green atom) from its neighbors which are inside the
cubic cell (shown in yellow) surrounding them. Here, the shaded red box is the actual MD
simulation box, and the white regions are image cells to simulate the bulk property.
2.6 Shifted potential

Figure 2.3: Actual and shifted interaction potential.
As discussed in the minimum image convention, atoms interact with their neighbors that
are enclosed within a cubic cell of dimension 𝑙 𝑥 ×𝑙 𝑦 ×𝑙 𝑧 . Because of this, interaction potential
could be discontinuous at the potential cut off distance. Hence, when the cell length is not large
enough then this truncation creates a significant discontinuity in the potential energy function.
16
Because of this discontinuity, we will not be able to differentiate the potential energy function, and
our derivation of forces will break down. To ensure that does not happened, each potential energy
function has a cutoff distance, rc, at which the value of the function becomes almost zero. Also, it
is important that the simulation box has to be greater than twice of this cut-off distance,
(max(𝐿 𝑥 ×𝐿 𝑦 ×𝐿 𝑧 )≥2𝑟 𝑐 ) , otherwise our force calculation will break down. Potential energy
function is also modified such at its value becomes zero at rc. The equation of the modified
potential energy function is shown in Equation 7, and Figure 3 shows the original and modified
potential energy curves. In general, rc is chosen in such a way that force contribution from atoms
at distance greater then rc is negligible. This choice of the value of the cutoff distance depends on
the material and functional form potential energy.
𝑢 ′
(𝑟 )={
𝑢 (𝑟 )−𝑢 (𝑟 𝑐 )−(𝑟 −𝑟 𝑐 )
𝑑𝑢
𝑑𝑟
|
𝑟 =𝑟 𝑐 𝑟 <𝑟 𝑐
0 𝑟 ≥𝑟 𝑐 Eq. 7
2.7 Parallel molecular dynamics
Using MD, we can simulate from million- to billion-atom systems. However, large scale molecular
dynamics simulations are intractable without a highly scalable parallel molecular dynamics (PMD)
framework. In PMD, the physical MD system is spatially decomposed into smaller domains of
equal volume and each of these domains are assigned to a separate processor. Figure 2.4(a-b)
shows the special decomposition of the MD simulation system into smaller cells and their
assignment to different processors. Here, processors are logically arranged according to the
topology of the physical sub-system, and atoms present in each sub-system are assigned to the
corresponding processors. Each processor is responsible for the computation of forces and
positions update of atoms assigned to it.
17

Figure 2.4: (a) Spatial domain decomposition of PMD simulation box. (b) Shows the assigned sub-
system to the processors along the three directions.
Spatial decomposition in PMD is done by giving each processor a unique index vector
𝑝 =(𝑝 𝑥 ,𝑝 𝑦 ,𝑝 𝑧 ) based on the Cartesian coordinate of the cell center in the physical system. These
index vectors are then mapped to a unique sequential id (𝑞 ) for each processor. Table 3 shows the
transformation between the sequential id and vector index of the processor. Here (𝑃 𝑥 ,𝑃 𝑦 ,𝑃 𝑧 ) are
total number of processers along x, y and z directions, respectively (Figure 2.4b).

Table 3: Algorithm to compute sequential and vector id of processors

The relation between sequential and vector id of processors
𝑝 𝑥 =
𝑞 𝑃 𝑦 ×𝑃 𝑧
𝑝 𝑦 =(
𝑞 𝑃 𝑧 )𝑚𝑜𝑑 𝑃 𝑦
𝑝 𝑧 =𝑝 𝑚𝑜𝑑 𝑃 𝑧
The calculation of sequential id from vector id
𝑞 =𝑝 𝑥 ×𝑃 𝑦 𝑃 𝑧 +𝑝 𝑦 ×𝑃 𝑧 +𝑝 𝑧
18
Since calculation of force on atoms requires their neighborhood information, we also need
to cache atoms from the neighbor processor (domain). Here, the thickness of the caching layer in
each direction has to be greater than the cut-off distance of the interaction potential. In two-
dimensional systems each processor has 8 neighbors as shown in Figure 2.5, and there are 26
neighbor domains in three-dimension. The most naïve method of atom caching will be the caching
of atoms of thickness Rcut from all 26 neighbors (face, edge and corner), which is called full shell
(FS) operation. However, using Newton’s third law, the number of neighbors required for caching
can be reduced to 13 in half-shell scheme (HS), and for the remaining 13 neighbors the computed
reactions will be send to them via Message Passing Interface (MPI) calls.

Figure 2.5: Here, each cell corresponds to a sub-system and are assigned to different processors
in parallel molecular dynamics (PMD). (a) Shows the catching layer regions (yellow) from the
neighboring processors for the highlighted red region. (b) Green arrows show migration of atoms
from one processor to another processor.
In PMD, resident atoms of each processors (sub-system) will move out of it and new atoms
will enter into it during each time step. Figure 2.5b shows a schematic of this atom migration
among processors. Hence at each time step, processers have to check for atoms that are migrated
out of the sub-system and migrated into it. For atoms that are migrated, processors have to remove
19
them from their assigned list and send the atoms to their new assigned processers via MPI calls.
For new atoms migrated into the processors, they will append the atoms into their atom list.
Since processors communicate amongst each other to keep track of migrated atoms. A
naïve implementation of commination pattern amongst processer will be where each processor
will send information to its neighbors and waits for their response before proceeding further. This
communication pattern will create a deadlock in certain situations, for example a cyclic loop where
the processors will wait for response from each other in an infinite loop (see Figure 2.6a). An
improved algorithm so as to avoid this deadlock situation will be a communication pattern which
consists of two phases: an even phase and an odd phase. During the even phase, processers with
even sequential id send information’s to the processers with odd sequential id, and during the odd
phase the processers with odd sequential id sends information to the processers with even
sequential id (see Figure 2.6b).

Figure 2.6: (a) Schematic showing communication deadlock. (b) Two-phase message passing
scheme to avoid deadlock.
20
2.8 Linked-list force computation

Figure 2.7: Schematic of linked list cell for force computation.
index 56 79 … 21 … 16 82 … 58 25
value 79 21 … 16 … 82 99 … 6 NULL
Table 4: Linked list data structure.
In brute force method, we compute forces on each atom by iterating over all the N atoms
present in a MD cell. This brute force method is an O(N
2
) algorithm and is computationally
inefficient for large values of N. A better choice for force computation is the linked-list based
algorithm that has computational complexity of O(N). In linked-list based method, physical sub-
system assigned to each processor, along with their cached layer is further sub-divided into smaller
bins and a linked-list data structure is used to organize the atom information present in each bin.
A schematic of the linked-list method is shown in Figure 2.7. This linked-list consist of a data
structure which contains the global id an atom and a pointer to the data-structure of next atom such
that all the atoms present in any specific bin can be traversed through this linked list (see table 4
and figure 2.7b). The algorithm consists of two steps: first, a linked list is created to traverse all
the atoms present in any specific bin. In the second step for each atom, we traverse through all the
21
atoms present in its resident bin and all its 13 neighbor bins and compute its interaction with those
atoms that are within the cutoff distance from its coordinate. Since each atom has approximately
constant amount of computation, k, from atomic pairs of its resident bin and neighbor bins the
computational complexity of this algorithm is O(kN). The reaction force on the neighbors from the
central atom are added to those atoms directly to them if they are not cached from the neighboring
processors otherwise its send to the remote machines where these atoms are physically present via
MPI calls.

22
Chapter 3: Mechanical Metamaterial
3.1 Background
Since prehistoric period, exploration and design of materials by mankind has given us
different classes of materials like metals, ceramics, polymers, semiconductors etc. As a result of
that, we are able to make advances in various fields of information storage, electronics, medical
devices and entertainment industry.
2,3
However, when we project all the materials discovered or
invented until now on a materials-property chart, we will observe that we have explored only a
tiny region of the entire space and still have many unexplored areas. Materials-property chart is a
multi-dimensional plot whose different axis represents properties of materials like strength,
stiffness, fracture toughness, thermal and electrical properties. Figure 3.1 shows the section of the
materials-property chart for strength vs. density for different classes of known materials. Material
scientist are continuously exploring techniques to design new materials so as to fill the gaps in the
materials-property chart. These techniques are in general based upon nanotechnology, chemistry
and microstructure of the materials. For example, chemistry-based approach involves developing
new alloys, polymers, and microstructure-based approach involves controlling the texture, defects
and distributing of different phases in the materials. Another approach to fill the gaps in materials-
property space includes creating hybrid materials, which are either composites or micro-
architecture materials. Composite materials are designed using combination of materials where
one phase has higher ductility and the other phase with higher strength whereas in micro-
architecture material, material and space is used to create hierarchical structures so as to design
new materials. From all these material design techniques, the micro-architectured materials have
been shown their novel properties that previously are not occupied in the materials-property
space.
28
Using these materials, we can design materials with properties like high-strength, high
23
stiffness and high fracture toughness at low density. These micro-architectured materials, which
are based on the scale of hierarchy and nodal connectivity can be tuned to provide materials with
unique combinations of strength, stiffness, and fracture toughness. This can also be used to design
materials with negative Poisson’s ratios and super-plasticity.

Figure 3.1: Shows the slice of the materials-property chart for all known materials.
In general, hierarchical solids represent materials that consist of internal sub-structures.
These hierarchy is usually defined by a number called hierarchical order number, which indicates
the levels of scale with identifiable structural pattern inside the material. The hierarchical order
number is 0 for usual solids and 2 for cellular structures such as honeycomb.
Hierarchy in structures is ubiquitous. Many natural and man-made materials exhibit such
structural hierarchy. A common example of hierarchical structure found in nature is bird nests. At
the macro level, bird nests are formed by leaves of grass, which are wavy and bend when nest
deforms. At the finer level, its constituents grass leaves are composed of three-dimensional form
like arrangement of cells with wavy sides. As a result of that, the cell walls of the grass leaf bends
24
under strain and make bird nest very damage tolerant under storms. Other examples of hierarchical
materials existing in nature are human bones and tendon, wood, and bamboo. For example, human
bone has rich hierarchical structures at different length scales. At a microstructural level, it has
hollow fibers composed of concentric lamellae and pores, which is known as osteons. These
lamellae are further built from fibers that comprise of sub-structure fibrils, which is further
composed from protein collagen. These hierarchy in structure the bone has several advantages in
terms of its property like for example, the composite structure of mineral microstructure controls
stiffness and protein fibers, osteons give slow creep deformation. The hierarchy of structural
design can be found even in some famous known monuments like the Eiffel tower, Garabit viaduct,
and Pompidou center.
Micro-architectured material
Micro-architectured or lattice material is a class of hybrid materials, which are made from
solid and space. Its solid part is made from struts and plates which are connected at the end of their
edges and give rise to an open structure. Lattice materials can be further divided into three groups
which are open-cell foam, closed-cell foam and lattice structures.
29
In open cell foams, the cells
within the structures are broken that allows the air to fill the spaces. They have very low density
and are sponge like structures. In closed cell foam, cells remain closed from the outside and hence
have no gaps for the air to fill the space inside the structure. Since no air is able to pass through
the foam that means they create an air tight barrier. The third-class lattice materials are lattice
structures which are similar to open cell foam but have higher nodal connectivity in comparison
of open and closed cell forms. Common examples of lattice structures are octet truss, kagome
lattice, and octrahedron.
25
The nodal connectivity inside lattice materials have important implication in their
deformation mechanism as low nodal connectivity favors bending-dominated deformation, and
higher nodal connectivity results into stretching-dominated deformation in the material.
30
The
rigidity of structures based on the numbers of struts and nodes it contains is determined by
Maxwell’s stability criterion.
29
It states that for a pin-joined frame containing b struts and j
frictionless joins to be both statically and kinematically determinate following equations has to be
satisfied for its rigidity which are 𝑀 =𝑏 −2𝑗 +3=0 in 2D and 𝑀 =𝑏 −3𝑗 +6=0 in 3D.
Here, statically determinate means that tension in all struts due to the external force can be
determined using the available number of independent equations and kinematically determinate
means that locations of all nodes can be uniquely determined using the length of each struts.
Structures that does not satisfy Maxwell Rule and have 𝑀 <0 are called mechanism. In this
system as shown in Figure 3.2a, frames have fewer beams than required for rigidity due to which
it has one or more degree of freedom that allow displacement along those directions. They behave
like mechanism for systems with pin-joins, and for rigid joined structures (open cell foam), they
show bending dominated deformation of its beams. Systems who satisfy the Maxwell’s criterion
show stretching dominated deformation for both pin-joined and rigid joined structures. Figure 3.2b
shows one such structure in which under external force the horizontal bar undergoes tension.
Systems with 𝑀 >0 have more beams than needed for its rigidity as shown in Figure 3.2c and
hence always remain in the state of self-stress. Self-stress means that beams remains in a state of
internal stress even in the absence of external force. In general, Maxwell Rule is a necessary
condition for rigidity as it does not take into consideration the mechanism and self-stress. A
generalization of Maxwell Rule in 3D is given by Calladine for the rigidity as 𝑀 =𝑏 −3𝑗 +6=
𝑠 −𝑚 . Here, s and m represent self-stress and mechanism present in the system. We can see from
26
this formulation that a system, Figure 3.2d, can have M equal to zero even when m and s are non-
zero.

Figure 3.2: 2D-frames showing Maxwell criteria (a) M<0, (b) M=0, (c) M>0 and (d) M=0.
In general, properties of lattice materials depend upon both their topology, relative density
and the properties of the base material from which it’s been constructed.
29
Figure 3.3 shows this
dependency of properties of lattice materials on various design attributes. Density of these
materials are dependent on the thickness of their struts which is quadratic for open cell foam, 𝜌 =
𝐴 (
𝑡 𝐿 )
2
and linear for closed cell form, 𝜌 =𝐴 (
𝑡 𝐿 ). Here, L is the length of the strut, t is its thickness
and A is constant of proportionality that depends on the topology of the lattice material.

Figure 3.3: Dependence of the property of the lattice materials on different parameters.

27
Mechanical Property of Micro-architectured material
Mechanical properties of lattice materials like Young’s modules E, yield strength 𝜎 and
Poisson ratio 𝜈 can be determined by simple beam theory in terms of their relative density (𝜌̃) and
the properties of their base material.
31,32
The Young’s modules, E, can be represented by a power-
law expression as
𝐸 𝐸 𝑠 =𝐵 𝜌̃
𝑏 . Here, 𝐸 𝑠 is the Young’s modules of the base material, and b is
constant of proportionality whose value is 1 for bending dominated structures and 1 for stretching
dominated structures. B depends on the geometry of the structure and is in between 0.3 to 0.5.
Yield strength (𝜎 ) of a defect free system can also be given in terms of the yield strength (𝜎 𝑌𝑆
) of
the base material as
𝜎 𝜎 𝑌𝑆
=𝐶 𝜌̃
𝑐 . In bending and stretching dominated material c is 1.5 and 1
respectively. The geometric factor C is again in between 0.3 to 0.5.
As mentioned before, open cell foams have low nodal connectivity and show bending
dominated deformation. In closed cell foams, cell walls are usually very thin, and they rupture or
buckle at very low stress. Because of this they contribute very little to the load bearing capacity of
the system and behave like an open cell foam. Hence closed cell forms also show bending-
dominated deformation. In comparison of these two systems, lattice structures (kagome, octet
truss) satisfy the Maxwell’s criteria and show stretching dominated deformation. Since stretching
dominated structures show linear dependence between 𝐸 /𝜎 and 𝜌̃, lattice materials have higher
strength and stiffness than foams for the same 𝜌̃. Also, due to this linear dependence between
strength and density, strength does not decrease rapidly with density, and we can synthesize
materials with very low density and high strength. Figure 3.4 shows the materials-property chart
for the modulus vs. relative density and strength vs. relative density for these two classes of
materials. We can see from Figure 3.4 that kagome and pyramidal lattices have higher modules
and yield-strength in comparison of foam.
28
In general, failure mechanism of bending-dominated structure is of three types, which are
plastic buckling, elastic buckling and crushing.
3,29
In elastic buckling, the ribs collapse is reversible
whereas in plastic buckling it’s irreversible and in crushing they fracture. Out of all these three
modes, the mode with minimum required load is usually the prevailing failure mechanism. Figure
3.5a shows the stress-strain curve for bending-dominated materials. It consists of a linear-elastic
regime up to its elastic limit, at which point the cell edges yields, buckle or fracture. The structure
continues to collapse at a nearly constant stress, 𝜎 , until the opposite sides of the cells impinge that
causes the densification of the system, and the stress goes up.

Figure 3.4: Comparison of modulus-density and strength-density chart in a stretching and bending
dominated structures.
In case of stretching-dominated structures, failure mechanisms are elastic-stretch
dominated, plastic-stretch dominated, buckling-dominated and stretch-fracture dominated.
Between plastic-stretch dominated and bucking-dominated deformation, bucking is the dominant
deformation mechanism for systems with slender struts as they usually buckle before they can
yield. Figure 3.5b shows the stress-strain curve for the stretching-dominated materials. In Figure
3.5b, we can observe a sharp decrease in stress after yielding. This happens due to the complex
deformation mechanism observed in these systems which are “hard” modes (tension, compression)
rather than “soft” ones (bending). Due to this reason the initial yield is followed by a plastic
29
buckling or brittle collapse of the ribs which leads to this post-yield softening in the stress-strain
curve. These stretching-dominated structures are in general good for light-weight structural
application due to their high strength and low-density. However, they are not suitable for energy-
absorbing applications, which requires stress-strain curves with a long flat plateau. Hence, the
major design challenge here is to reduce this post-yield softening of these structures.

Figure 3.5: Stress-strain curve for (a) a bending dominated and (b) a stretching dominated
materials.
A major difference between two known lattice structures - octet truss and kagome lattice -
is that octet truss is stretching dominated and kagome lattice is stretching periodic-bending
dominated structure. For example, in 2D kagome lattices, stacked triangular unit cells are
stretching-dominated however due to the tessellation pattern of the structure, the resulting lattice
shows bending-dominated deformation at the nodes. Because of this unusual deformation
mechanism, their fracture toughness scales as 𝐾 𝐼𝐶
∝𝜎 𝑌 𝜌̃
0.5
whereas for octet truss it scales linearly
with 𝜌̃. 3D kagome lattice can exploit this deformation mechanism so as to create materials with
enhanced fracture toughness however their fabrication is still a challenge.
33

30

Fabrication of Micro-architectured materials
Techniques known for the fabrication of these micro-architectured materials are limited and are in
general based on photolithography, 3D printing, self-propagating polymer waveguide and two-
photon direct laser writing lithography.
34-36
In photolithography, 3D lattice materials are created
by first creating a 2D-lattice using soft-lithography and electrodeposition, which is then folded to
create 3D structure. Hence, it is not suitable to create very complex and intricate 3D-truss
structures. 3D printing (3DP) is another approach as it allows us the complete control over the
geometry of the lattice structures. In 3DP, structures are formed by rastering process, which
deposits adhesion/bulk materials to create thin cross-sectional areas. The process starts from the
bottom and in incremental fashion cross-sectional areas are built until the structure is complete.
However, the major drawback of this method is that it is very time consuming, and the resolution
of the created structures are in the order of 100 𝜇𝑚 . Hence, the most commonly used technique
for lattice material fabrication are based on self-propagating polymer waveguide and two-photon
direct laser writing lithography as they allow us to create micro-lattices in the range of 𝜇𝑚 to 𝑛𝑚

Figure 3.6: Fabrication of micro-lattice by self-propagating polymer waveguide.
4

Self-propagating polymer waveguide utilizes the self-propagating phenomena of
photosensitive monomers to create 3-dimensional truss structurers. The schematic of this process
is shown in Figure 3.6. Figure 3.6a shows the creation of the polymer templet of the truss structure.
31
Single point exposure of light in the photosensitive monomer creates self-propagating polymer
waveguide inside it. This happens due to the tunneling effect of the incident light inside the
monomers as the waveguides are formed. In general, collated ultraviolet (UV) rays are used as an
incident beam, and they are passed through a patterned mask into the monomers. When UV rays
propagates through the monomers via these patterned masks, they create multiple waveguide
inside the monomer and forms the a 3-dimensional polymer templet of the truss. The uncured
monomers are removed after this step. To create a metallic micro-lattice, metallic film is deposited
on these polymer templets via electroless plating as shown in Figure 3.6b. In the last step, polymer
templet is etched out from the structure, and we get a micro-lattice where each strut are made up
of hollow tubes (Figure 3.6c). During this process, the thickness of the hollow tubes in the micro-
lattice is controlled by the electroless plating’s time. Schaelder et al. has created an ultralight
cellular lattice of kagome made from polycrystalline nickel, where each strut had diameter (100-
500 𝜇𝑚 ) and wall thickness (100-500 𝑛𝑚 ).
4
They have used thiol-ene liquid monomer to create
the micro-lattice templet, and the conformal nickel-phosphorus thin film is deposed on it using
electroless plating. For electroless plating, they have used nickel sulfate as nickel source, sodium
hypophosphite as reducing agent and sodium malate and acetic acid as complexing again. The pH
of the bath is kept at 4.9 by the addition of ammonium hydroxide and the whole process is
performed at 80°C. 3M NaOH is used an etching agent, which is performed at 60°C for 20 hours.
They showed that these micro-lattices have very low density (~0.9mg/cc
3
) and high strength, and
they can be compressed up to 50% with complete recovery. Their deformation mechanism is
bending-dominated, which happens due to the absence of struts in the basal plane of these
structures and smaller wall thickness of struts as both of them favors bending dominated
deformation.
32

Figure 3.7: Fabrication of micro-lattice by two-photon direct laser writing lithography.
34

The last technique, which is two-photon direct laser writing lithography (TLP) can create
lattice structures with complex shape for length scale between 𝜇𝑚 to 𝑛𝑚 at very high precision.
Montemayor et al. have used this technique to create micro-lattices of octet truss, kagome and
octahedron, where the base material was either alumina or gold. A schematic of TLP process is
shown in Figure 3.7. At first a CAD design of the micro-lattice is created and after that its vertices
are transferred to the NanoWrite Software in the form of set of points. NanoWrite is an interface
with communicates with the Nanoscribe TPL instrument who creates the polymer templet. After
that, IP-Dip photoresist which is specially designed by Nanoscribe GmbH to optimize the TPL
process is exposed to 780 nm femtosecond pulsed laser. These lasers are focused to a region, which
is called voxel that are elliptical in shape and have sufficient energy to initiate cross-linking in the
photoresist material to create the polymer templet. These voxels are traced according to the
geometry information provided by the CAD design and creates a 3-dimensional polymer templet.
The laser power and writing speed is very critical for this process and are in general optimized for
33
different structures and materials. For example, with low writing speed and laser power the
structures collapse due to the capillary force whereas with high write speed the structure is
destroyed. The created polymer templet is developed for 30 mins in propylene glycol mono-methyl
ether acetate (PGMEA) before metallic coating. Metallic coating is done on these templets using
sputter deposition and after that these polymer templets are etched out from the micro-lattice by
O2 plasma treatment. Before the O2 plasma treatment, the edges of the nanolattice are milled on
two sides using focus ion beam so as to expose the polymers for O2 plasma treatment. The resulting
micro-lattice after this step are made up of hollow struts of the base materials. To test the thickness
of the of hollow struts, the fabricated structures are sliced 1/3, 1/2 or 2/3 along the depth of the
lattice using focused ion beam and inspected visually. Meza et al. have created octet truss of
amorphous alumina by this method with densities spanning between 6.1 to 249 km/m
3
. Their
density was much smaller than that of bulk alumina (~2900 kg/m
3
).
6,7
The created structures
showed density dependent deformation mechanism, where thick walled lattice with length to
diameter ratio greater than 0.3 collapsed due to brittle fracture of the tube walls. Whereas lattices
with thin walled alumina tubes deformed via shell buckling. These materials were compressed up
to 50% strain with complete recovery and have elastic failure.
Current Work
Lattice materials have tunable property due to which they can be used for energy
application, structural design, amours, thermal insulator, railways, etc. While strength is the
materials property which can we tuned by this hierarchical structural design, other bulk properties
like toughness is highly depended on the weakest path in the system, for example, crack length.
Hence, designing materials with high-toughness, high-strength and low density is still a challenge.
34
Another important aspect of these hierarchical material which has not been extensively studied is
the deformation behavior of hollow lattices at atomic level, and how is their atomistic deformation
mechanism is different from the lattice structure where constituent strut are nanorods instead of
hollow tubes. In this work, we have studied the deformation mechanism of solid and hollow
kagome lattices made from nickel. We have also studied the mechanism of mechanical collapse of
their constituent struts, which are nickel nanorods and hollow tubes.
3.2 Simulation details
3.2.1 EAM interaction potential
The Embedded Atom model (EAM) is an empirical force filed, which is suitable to model
metallic systems, in particular FCC crystals like nickel, copper, etc.
37,38
It is based on the concept
of local electron density and is able to describe bonding in metallic systems as it can account for
the dependence of the strength of the individual bonds on the local environment. Hence, it can
describe surfaces, defects and plastic deformation by dislocation motion in metallic systems.
Another advantage of EAM potential over pair interaction potential is that it can model variation
in bond strength with coordination. It means that with increase in coordination decreases the
strength of each of the individual bonds, which increases the bond length.
In addition to the EAM potential, many different methods have been developed to describe
the metallic systems, for example, Finnis-Sinclair potential, glue model, effective medium theory,
and corrected effective medium theory (CEM). All these models are based on physical arguments.
For example, tight-binding model, effective medium theory, etc. But all of them have similar
expression of the potential energy. The functional form of the EAM potential is given in equation
3.1 for a system consisting of N atoms. It consists of two terms, a pair-wise potential function, ij,
between atoms and an embedding function F, which represents the energy required to place an
35
atom into an electron cloud. This embedding function depends on the local electron density around
an atom.
𝑈 =∑ (𝐹 𝑖 (𝜌 𝑖 )+
1
2
∑ 𝜙 𝑖𝑗
(𝑟 𝑖𝑗
)
𝑖 ≠𝑗 )
𝑁 𝑖 Eq3. 1
Where
𝜌 𝑖 =∑ 𝑓 (𝑟 𝑖𝑗
)
𝑗 Eq3.2
The exact functional form and interpretation of F, f,  depends on the method under
consideration. From the point of view of effective medium theory or the embedded atom method,
the energy of an atom is determined by the local electron density at the position of the atom, where
f describes the contribution to the electron density at the site of the atom due to all its neighbors.
Hence the summation of the function f around an atom is a measure of local electron density around
that atom. As mentioned above, embedding energy F is the energy associated with placing an atom
in an electron environment described by 𝜌 and  energy associated with pair interaction. In
equation 3.1-3.2, rij is the distance between central atom i and its neighbor j. The parameters of the
EAM potential is usually fitted using the DFT results and experimental data. Beside defects and
surface properties, EAM potentials is able to give correct value for of the lattice constant, cohesive
energy, bulk modulus, elastic constant and melting point. For EAM potential, the force on atom is
given by equation 3.3-3.4.
𝐹 𝑖 ⃗⃗
= −Δ
⃗⃗
𝑟 𝑖 ⃗⃗⃗
𝑈 =−Δ
⃗⃗
𝑟 𝑖 ⃗⃗⃗
[∑ (𝐹 𝑖 (𝜌 𝑖 )+
1
2
∑ 𝜙 𝑖𝑗
(𝑟 𝑖𝑗
)
𝑖 ≠𝑗 )
𝑁 𝑖 ] Eq3. 3
𝐹 𝑖 ⃗⃗
=−∑ [
𝜕 𝐹 𝑖 (𝜌 )
𝜕𝜌
|
𝜌 =𝜌 𝑖 ×
𝜕𝑓 (𝑟 )
𝜕𝑟
|
𝑟 =𝑟 𝑖𝑗
+
𝜕 𝐹 𝑖 (𝜌 )
𝜕𝜌
|
𝜌 =𝜌 𝑗 ×
𝜕𝑓 (𝑟 )
𝜕𝑟
|
𝑟 =𝑟 𝑖𝑗
+
𝜕 𝜙 𝑖𝑗
(𝑟 )
𝜕𝑟
|
𝑟 =𝑟 𝑖𝑗
]
𝑖 ≠𝑗
𝑟 𝑖 ⃗⃗⃗ −𝑟 𝑗 ⃗⃗⃗
𝑟 𝑖𝑗
Eq3. 4
In this work, to simulate the flat punch compressive of Ni nanorod, nanotube and kagome
36
lattice we have used EAM potential of Ni developed by Mishin for NiAl system.
39
The properties
Ni predicted by EAM potential along with the DFT/ experimental data is given in Table 3.

Properties DFT/Experimental value EAM potential
Lattice consent 3.52 A 3.52 A
Cohesive energy -4.45 ev/atom -4.45 ev/atom
Bulk modulus, B 1.81 Pa 1.81 Pa
Elastic constant, C11 2.47 Pa 2.47 Pa
Elastic constant, C12 1.47 Pa 1.48 Pa
Elastic constant, C44 1.25 Pa 1.25 Pa
Surface energy (110) 2280 mJ/m
2
2049 mJ/m
2

Surface energy (100) 2280 mJ/m
2
1878 mJ/m
2

Surface energy (111) 2280 mJ/m
2
1629 mJ/m
2

Defect, Stacking fault 125 mJ/m
2
125 mJ/m
2

Defect, grain boundary (210) 1572 mJ/m
2

Defect, grain boundary (310) 1469 mJ/m
2

Vacancy, 𝐸 𝑣 𝑓
1.60 eV 1.60 eV
Interstitial, 𝐸 𝐼 𝑓
5.86 eV
Table 3.1: Comparison of physical and mechanical properties of Ni predicted by EAM potential
with the experimental data
3.2.2 Defect identification in FCC crystal
Mechanical deformation of material by molecular dynamics simulation forms complex
structures inside the materials. For examples, surfaces, point defect, mobile and immobile
37
dislocations, grain boundaries and twin boundaries. Simple coordination number-based analysis
of atoms will not be able to identify these structures as these are complex functions of local
geometry around the defect atoms. Three most commonly used techniques to identify dislocations
in monoatomic FCC crystals are centrosymmetry parameter (CSP), common neighbor analysis
(CNA) and common neighborhood parameter (CNP).
40

Centrosymmetry parameter (CSP)
The CSP method is developed by Kelchner, Plimpton and Hamilton and it based on the
fact that centrosymmetric material remains centrosymmetric under homogeneous elastic
deformation. Hence in such materials, nearest neighbors of an atoms form certain pairs of equal
and opposite bonds around it. In case of FCC crystals, each atom has 12 nearest neighbors and
they forms 6 pairs of equal and opposite bonds around it. Hence, the CSP parameter for an atom
in FCC crystal is zero and any atoms whose value is greater than zero is considered as defect. The
equation to compute CSP parameter is given in equation 3.5 where the sum is over the opposite
pairs of neighbors around the central atom, i. A schematic of equation 3.5 is given in Figure 3.8b.
𝐶𝑆𝑃 =∑ |𝑅 𝑖 ⃗⃗⃗
+𝑅 𝑖 +6
⃗⃗⃗⃗⃗⃗⃗⃗
|
2
6
𝑖 =1
Eq3. 5
Common Neighbor Analysis (CNA)
Common Neighbor analysis (CNA) is developed by Honetcutt and Andersen in which atoms
are descried by a set of four parameters described below.
i) The first parameter has value 1 or 2 depending up on whether the two atoms, 𝛼 and 𝛽 ,
are each other nearest neighbors (1 if they are neighbor and 2 otherwise).
ii) It represents the number of common nearest neighbors, k, 𝛼 and 𝛽 have.
iii) The third parameter indicates the total number of bonding between the common
neighbors of 𝛼 and 𝛽 .
38
iv) Distinguishes atom from each other when they have same values for i-iii but have
different order of bonding among its nearest neighbors.
In FCC crystal, CNP parameter for an atom is (1421) which means each atom and its nearest
neighbor has four common neighbors amongst each other and between these common neighbors
we have two bonds. A schematic of CNP parameter in FCC crystal is shown in Figure 3.8a.

Figure 3.8: Schematic of CNA, CSP and CNP parameters for the central brown atom. Here, CNA
of the central atom in (a) is (1421), (b) opposite neighbors of the central atom are shown in same
color in for CSP and (c) common neighbor of the neighbors of the central atom is shown in yellow
for CNP.
40

Common neighborhood parameter (CNP)
Common neighborhood parameter (CNP) is developed by H. Tsuzuki, P. S. Branicio and
J. Rino. CNA gives a single number for each atom, which can used to distinguish different types
of defects in FCC, HCP and BCC crystal. CNP is based on the concept of CSP and CNA. CNP
value is determined by the total number of nearest neighbors of an atom and common neighbors
between the atom and its nearest neighbors.
The functional form of CNP (Qi) is shown in equation 3.6. In equation 3.6, nj is the total
number of nearest neighbors of atom i and njk is the total number of nearest neighbor common
between atom i and j. Rij and Rik is the distance between the central atom i and its nearest neighbor
39
j with their common neighbor k. A schematic of the CNP parameter is shown in figure 3.8c and
its value in FCC crystal structure is 0.
𝑄 𝑖 =
1
𝑛 𝑖 ∑ |∑ (𝑅 𝑗𝑘
+𝑅 𝑖𝑘
)
𝑛 𝑖𝑗 𝑘 =1
|
2
𝑛 𝑖 𝑗 =1
Eq3. 6
All these CSP, CNA and CNP parameters can be used identity defects and dislocations in
other monoatomic systems like BCC and HCP crystal structures. A comparison of these
parameters in different crystal structure is given in Table 3.2. For diatomic system like NiAl, these
parameters in their native form does not work, however variation of CNP parameter exists to
identify defects in the dia-atomic systems. Also, these parameters can only identify whether an
atom is a defect or a dislocation, however it does not give us any information about the burgers
vector of the dislocation.
Crystal Structure CSP CNA CNP
FCC 0 A
2
0 A
2
1421
BCC 0 A
2
0 A
2
1421, 1422
HCP 6.6 A
2
4.4 A
2
1441, 1661
Table 3.2: Comparison of CSP, CAN and CNP values for FCC, BCC and HCP crystal structure.
Dislocation extraction algorithm (DXA)
Dislocation extraction algorithm (DXA) is developed by A. Stukowski and K. Albe, which
is based on CNA.
41
It can identify different types of dislocations present in FCC crystal along with
burgers vector, stacking fault and twin boundaries. It consists of three steps. First, it calculates the
CAN parameter for each atom so as to distinguish between the crystalline atoms and defects. In
step two, it constructs a closed, orientable, two-dimensional manifold (an interface mesh) that
separates crystalline atoms from defects. In step three, for each dislocation segment, it generates a
40
Burgers circuit on this manifold and advances it both directions to captures the dislocation shape.
The output of the DXA algorithm is dislocation lines and its burgers vector.
3.3 Uniaxial compression of nickel nanorod and nanotube
3.3.1 Simulation setup
Uniaxial compression of nickel nanorod and nanotube is studied using homogeneous and
flat punch compression for two different aspect ratios, which are 1:4 and 1:8.
42-47
The total number
of atoms in these systems are between 270,000 to 2,500,000. Ni nanorod of diameter 10 nm is
created from bulk crystalline nickel such that its length to diameter ratio is either 1:4 or 1:8. After
that nanotubes of thickness 2 nm is created by removing inner cylindrical region of diameter 8nm
from the nanorods. In case of flat punch compression, additional circular disk of diameter 12nm
and thickness 5 nm is placed at the top and bottom surfaces of the nanotubes and the nanorods.
Here, the bottom circular disk acts as the substrate which is kept fixed during the compression and
the top circular disk is moved at constant speed and it’s called a flat punch. Schematic
representation of homogeneous and flat punch compression is shown in Figure 3.9.

Figure 3.9: Schematic representation of homogenous and flat punch compression.
After creating the systems, they are relaxed by setting the atomic velocities zero in every
10 steps for 10,000 steps. The time step used here is 2.5 fs. These systems are further relaxed for
50 ps using steepest descent method. After relaxing the systems, they are heated slowly in steps of
41
50,000 until the system temperature reaches 300 K, then thermalized at 300K for 400ps under
NVE ensemble. The relaxed configurations of nanorod and nanotube are subjected to compressive
load via homogeneous and flat punch compression. In case of homogeneous compression, the
system is uniformly compressed by 1% in 4,000 steps by rescaling the atomic coordinates along
[001] direction. Whereas in case of flat punch compression, the flat punch is moved along negative
[001] direction by 4 Å in 4,000 steps and the substrate is kept fixed during the compressive loading.
Here the speed of the flat punch is 40 m/sec and the strain rate applied on the system in both cases
(homogeneous and flat punch) is approximately 2×10
8
sec
-1
. After compressing the systems, they
are relaxed for 400 ps under NVE ensemble. In case of flat punch, both the substrate and the flat
punch atoms is kept fixed during the relaxation phase. The set of compression and relaxation steps
is repeated until the mechanical collapse of the systems. Periodic boundary condition is applied
only along (001) direction in case of homogeneous compression which is the direction along which
the systems are compressed.
3.3.2 Stress-strain curve
The stress–strain curve for homogenous and flat punch compression of Ni nanorod and
nanotube with aspect ratio 1:4 and 1:8 is shown in Figure 3.10. We can observe from the stress-
strain curve that the onset of plastic deformation in nanorod and nanotube occurs at around 7%
compression during homogenous compression. However, in case of the flat punch compression
mechanical collapse happens at around 3.8% compression and 3.3 strains respectively for systems
with aspect ratio 1:4 and 1:8. Also, the mechanical collapse of nanotube is observed at slightly
lower strain than in the case of the nanorod system.
The reason for the mechanical collapse of the systems at higher strain in case of
homogenous deformation can be explained in terms of homogeneous and heterogeneous
42
nucleation of defects in these systems. In case of flat punch deformation, the contact area of the
nanotube/nanorod with the substrate or the flat punch act as a heterogeneous nucleation site for
dislocation since the nanotube/nanorod experiences forces both from the substrate and the flat
punch. As a result, the stress concentration is higher in those regions. However, in case of
homogenous compression the deformation happens at higher strain due to the absence of such
dislocation nucleation site. At higher strain, the critical resolved shear stress for dislocations to
nucleate exceeds at multiple locations inside the system, which results into the homogenous
nucleation of dislocation at multiple locations, and subsequent plastic deformation inside the
system.

Figure 3.10: Stress-Strain Curve for homogeneous and flat punch deformation of nanorod and
nanotube of aspect ratio 1:4.
To understand the deformation mechanism in these system during their mechanical
collaspe, we have analysed the stress distribution, nucleation and propagation of dislocation inside
the during their failure. In all systems during flat punch compression, it is observed that
43
deformation happens due to the nucleation of 1/6<112> Shockley partial dislocation near the
substrate or flat punch. However, the motion of dislocations and their interaction amongst each
other is very different in nanorod and nanotube.
3.3.3 Mechanical collapse of nickel nanorod

Figure 3.11: Systematic representation of dislocation motion in nanorod. Shockley partial
dislocations are shown in yellow and hirth lock is shown in red. Various {111} planes on which
dislocation motion is happening are shown in blue, orange and green color.
Figure 3.11 shows the nucleation of multiple Shockley partial dislocation on {111} plane
from the surface of nanorod at the contact region between nanorod and the substrate during the
onset of plastic deformation in nanorod with aspect ratios 1:4. As these dislocations start
propagating on different (111) planes, they start interacting with other partial dislocations. The
interactions of these Shockley partial dislocations lead to the formation of immobile strain rod
dislocations with burgers vector 1/3<010> by the reactions as shown in equation 3.7. These strain
rod dislocations form an immobile hirth lock as shown in figure 3.11b in red color.
48
Due to this,
44
further dislocation motion advances in only two of the {111} plane (blue and green planes in
Figure 3.11c) as dislocation glide on other planes get blocked.
1
3
[010]=
1
6
[112]+
1
6
[1
̅
12
̅
] Eq3. 7
Eventually the partial dislocation on these plane reaches the surface of the nanorod and
gets annihilated. This motion of single dislocation through the nanorod creates a stacking fault in
the system (See Figure3.15). Further deformation of the system happens locally by dislocation
nucleation only on one of the (111) planes (shown in blue in figure 3.11c), since the path of the
dislocation glide on all other (111) planes are blocked by this one. We can observe from Figure
3.11d that at the interaction of different glide planes, dislocation reactions lead to formation of
hirth lock (red line at the intersection of blue and green planes). Along with the surface of nanorod,
these hirth lock also acts as the nucleation site for dislocation. Similar dislocation nucleation
mechanism is observed in materials by single-arm nucleation. We can see from Figure 3.11d-e that
partial dislocation loops (yellow) are nucleating both from the surface and also from the hirth locks
(red). Eventually, the hirth lock breaks down into Shockley partial dislocation again and their
motion leads to the removal of other two stacking fault (green and the orange one) from the system
as shown in Figure 3.11f. After that the nucleation and propagation of dislocations continues on
only two parallel {111} surface (blue in Figure 3.11f) where on each surface one end acts as a
source and the other end acts as sink for the dislocations.
This mechanism of deformation in nanorod can be explained by dislocation starvation
mechanism proposed earlier as deformation mechanism for nano-scale materials.
49
Unlike micron
size sample, in nanorods dislocations have very small distance to move before they get annihilated
at the surface. Because of these they never get chance to multiply and hence we have very low
dislocation density in the system. System always remains in a dislocation starved state and
45
deformation of the nanorod happens only due to these sudden bursts of dislocations at the nanorod
surface. They glide though the nanorod and leave the nanorod at the opposite surface. As
dislocations nucleate and propagate through the nanorod, the overall stress in the system decreases
and then stops further nucleation of dislocations. Because of the stress release event, we observe
the drop of stress level in the systems during its mechanical collapse. Figure 3.14 shows the
dislocation density and stress in the system as a function of time during homogeneous compression
using this the relation between dislocation nucleation and stress distribution in the system is
manifested. Here dislocation density and the burgers vector of the dislocation has been calculated
using dislocation extraction algorithm.
3.3.4 Mechanical collapse of nickel nanotube

Figure 3.12: Systematic representation of dislocation motion nanotube. Shockley partial
dislocations are shown in red, yellow, green and magenta.
In case of nanotube, deformation happens again due to the sudden burst of dislocations in
the nanotube and as the dislocation activity continues in the systems, stress in the system drops
and finally deformation of the nanotube stops. However, the major difference here is in the
mechanism of dislocations nucleate and propagate inside the nanotube. In nanotube, as shown in
Figure 3.12a, a single 1/6<112> Shockley partial dislocation nucleates from the outer wall. After
46
that this dislocation glide on the (111) and reaches the inner wall. Here, only a part of the
dislocation gets annihilated at the inner wall and its remaining segment split into two separate
dislocations as shown in yellow and green in Figure 3.12b. These two dislocations then move
separately in the nanotube as shown in Figure 3.12b and causes further increase in length of the
stacking fault inside the nanotube, which is created inside the nanotube due to their motion. Finally,
they meet again and annihilate each other as shown in Figure 3.12c and create a complete stacking
fault in the system.
When we look at the trace of the motion of this dislocation at the onset of mechanical
collapse in Figure 3.11a-c for nanorod and compare that with of nanotube in Figure 3.12a-c, we
observe that dislocation motion follows a straight line in case of nanorod from its nucleation to
annihilation site. However, in case of nanotube it follows a curved trajectory. Further deformation
of the nanorod happens due to the spiral motion of new dislocations nucleated on this path as
shown in Figure 3.12 d-f, first due to sudden burst of partial dislocation and then due to series of
single dislocations. Figure 3.12d shows the sudden outburst of dislocation nucleation at the inner
and outer wall of the nanotube at multiple places in the system. Because of this sudden burst of
multiple of dislocations nucleation what happens is that after a dislocation split into two
independent segments, only one of its segments get annihilated immediately (Figure 3.12d-f). For
example, we can see in Figure 3.12d that a dislocation nucleated from outer surface (shown in
yellow) split into two independent segment and at the same time another dislocation nucleate
nearby (shown in magenta in figure 3.12e). Figure 3.12f shows the annihilation of one of
independent dislocation segments of these yellow and magenta dislocations by each other (shown
in green). Thus, we observe here that lifetime of some dislocation segments is smaller compared
to that of others in nanotube. Due to which they effectively cause very little deformation of the
47
system. Now, the deformation of the nanotube happens only due to the surviving segments of the
dislocations and they follow a curved path (Figure 3.12f). Since this dislocation has to follow a
curved path hence they need a longer time to create a complete slip step in the system. Also, the
nanotube does not have any specific site which acts as a sink for these dislocations. Dislocations
can annihilate only when one dislocation meets with another dislocation whose direction of motion
is opposite to that of its direction. They merge and leave the nanotube either from inner or outer
surface. Due to the anisotropy in the motion of total number of dislocations moving in opposite
direction, lifetime of dislocations in nanotube is larger compared to that of nanorod as shown in
Figure 3.14, and they also require longer time to create a single slip in the system. The deformation
mechanism in nanorod and nanotube with aspect ratio 1:8 is exactly same as that of above systems.

Figure 3.13: Systematic representation of dislocation motion in nanorod (a-c) and in nanotube (d-
f) during compression at mechanical collapse.
The mechanism of deformation in nanorod and nanotube in case of homogenous
compression is same as that of flat punch compression however, here dislocation nucleation
happens at multiple sites throughout the system. Due to this deformation mechanism, both nanorod
and nanotube shows severe plastic deformation after the mechanical collapse at multiple locations
48
which is not the case in flat punch compression where the deformation happens locally near the
substrate. The mechanism of deformation of the nanorod and nanotube under homogeneous
compression at the onset of mechanical collapse is shown in Figure 3.13. It can be seen from Figure
3.13 that multiple dislocations nucleate simultaneously at different locations on the surface of
nanorod (Figure 3.13a) and nanotube (Figure 3.13d). Here, we also observe that even though we
have multiple deformation site in the systems, in nanorod dislocation lines follows a straight-line
path and in nanotube it follows a curved path similar to that of flat punch compression. Due to the
deformation of the system at multiple sites, both nanorod and nanotube are heavily deformed and
shows multiple twin regions after the collapse (Figure 3.13 c and f). This is not the case with flat
punch compression where damage happens locally in the system near the substrate.
3.3.5 Comparison of mechanical collapse in nanorod vs. nanotube
Figure 3.14 shows the dislocation density in nanorod and nanotube of aspect ratio 1:4 in
case of homogeneous compression during mechanical collapse. We can see from the plot that
decay of the dislocation density in nanotube is slower than of nanorod. This happens because
dislocations have to travel larger distance to create a complete slip step in the nanorod and
dislocation annihilation happens here only when two dislocations moving in opposite directions
meet each other and forms a loop which then annihilate at the surface and hence takes longer time.

49

Figure 3.14: Dislocation density and stress drop in the system as a function of time during the
mechanical collapse of nanorod and nanotube in homogeneous compression.

Figure 3.15: Formation of twin region inside the deformed regions of nanotube and nanorod due
to the dislocation motion.
Another scientific question would be what happens locally to the structure of nanorod and
nanotube in the deformed region of these system and which of these systems has deformed the
most. To obtain atomistic level insight, we have analyzed the structure of the deformed region in
these systems. In summary, the propagation of a single Shockley partial dislocation leads to the
formation of stacking fault in the systems (Figure 3.15). The deformation of FCC metals via
twinning or slip is reported as deformation mechanism in nanorods under uniaxial tension and
50
compression.
50,51
Here, both in nanorod and nanotube due to the repeated nucleation and
propagation of <112> partial dislocations at the exact same region at adjacent parallel {111} planes
(blue planes in Figure 3.15) leads to the formation of twin region in the system as shown in Figure
3.15b. The width of these twin increases further with dislocation activity in that region. Hence to
quantify the extent of damage in nanorod and nanotube, the width of this twin region is a good
candidate. We have calculated the width of the twin region in all these four systems after the
mechanical collapse of the systems. The width of the twin region is taken as the perpendicular
distance between its boundary as shown in Figure 3.15c in red. The width of this twin region for
nanorod with aspect ratio 1:4 and 1:8 is 27Å and 31Å respectively. While for nanotube its 16.5Å
(aspect ratio 1:4) and 28.6Å (aspect ratio 1:8). We can see that width of the twin region in nanotube
is smaller compared to that of nanorod for both systems. This can again be explained in terms of
dislocation motion in these systems. Firstly, to increase the width of the twin region, a dislocation
line has to glide through the entire cross-section of the nanorod or nanotube. As we have seen
earlier that in nanorod the path of the dislocation lines is straight. Also, the presence of hirth lock
which acts as additional nucleation site which further decreases the time required to increase the
width of the twin region. Since dislocation nucleates independently on both side of the hirth lock
and they only have to travel half the distance to increase the width of the twin region. We also
observe here that dislocation lines are curved and larger in nanorod hence they swipe larger area
of the nanorod in each time step. In case of nanotube, due to the curved path of motion of the
dislocations, it takes larger time to increase the width of the twin region. Also, here lots of
dislocations lines get annihilated very quickly and they result in effectively no deformation of the
system. To validate this hypothesis, we have calculated the time required for first few dislocations
51
to create a complete slip in the system at the onset of mechanical as shown in table 3.3. We can
see that in nanotubes dislocations takes almost double time to increase the width of the twin region.
Dislocation 1 Dislocation 2 Dislocation 3
Nanotube 18.75ps 16.25ps 15.25ps
Nanorod 9ps 8.5ps 9.5ps
Table 3.3: Time taken by three dislocations to create a slip step
3.3.6 Summary
In summary, mechanical collapse in homogenous compression happens at higher strain
then that of flat punch compression. In homogenous compression, mechanical collapse of nanorod
and nanotube happens at multiple location due to homogenous nucleation of dislocation
throughout the sample. In flat punch compression, amount of damage is smaller compared to that
of homogenous deformation and it happens near the substrate. During flat punch compression, the
extent of damage in nanotube is smaller compared to that of nanorod. This happens due of fact that
dislocation moves along a curved path in nanotube and hence requires larger time to increase the
width of the twin region in the system. While in nanorod dislocation segments are larger in length
and moves along a straight line and hence requires less time to increase the width of the twin
region.
3.4 Flat punch compression of nickel kagomé lattice
3.4.1 System setup
The deformation mechanism of the kagomé unit cells (KUCs) at the atomistic level is not
well understood, especially the response of these system as a function of their density and wall
thickness of their struts. To understand the deformation mechanism of these system at the atomic
level, we have performed flat punch compression of KUCs made from Ni nanorods and hollow
52
tubes by molecular dynamics (MD) simulation. To perform the MD simulations, the KUC of
dimension 32.5𝑛𝑚 × 32.5𝑛𝑚 × 32 .5𝑛𝑚 is created from Ni cubic lattice of by cutting eight
cylinders from it such that their axis is aligned along the four diagonals of the cubic lattice. These
cylinders have the outer diameter 10 nm and are either hollow (with thickness 2nm) or solid,
depending upon whether the constituting element of the KUC is nanorod or hollow tubes. The
length the struts in each KUC is 25 nm, and the angle between them at the node is 45
0
. Hence, the
overall density for the systems made from nanorod and hollow tubes is 2.67 gm/cc and 4.34 gm/cc,
respectively. To perform uniaxial compression, a flat punch and a fixed substrate of dimension
52.8𝑛𝑚 × 52.8𝑛𝑚 × 2 𝑛𝑚 is also added at the top and bottom of the KUCs such that initial
gap between them is 3Å. Figure 3.16 shows a schematic of the KUC, which shows its 8 struts and
the nodal region.
The entire system is first relaxed using conjugate gradient, and then heated to 25 °C using
the NVT ensemble in 50,000 MD steps. Then it is relaxed at 25 °C for another 100 ps in the NVE
ensemble. This process creates the relaxed configuration for the flat punch compression
simulations. The compression of the KUCs consists of two phases, a compression phase, which is
followed by a relaxation phase. During the compression phase, the flat punch is moved by 1.5 Å
in 30,000 steps along the negative y-axis while the substrate is kept fixed (see figure 3.16). In the
later phase, the kagomé unit cell is relaxed for 80,000 steps under the NVE ensemble, and during
this time both the substrate and the flat punch is kept fixed. The entire process is repeated again
until, the KUC is compressed up to 20%. This completes the loading phase of the compression,
which is followed by the unloading phase. During the unloading phase, the direction of motion of
the flat punch is switched to positive y-axis, and the entire cycle is repeated again until, the stress
53
inside the KUC drops to zero. The strain rate used here during the loading and the unloading phase
is 1.6×10
−7
sec
−1
, and 2.5 fs as the time step for simulation.

Figure 3.16: Schematic of flat punch compression of KUC.
To understand the deformation mechanism in structures made from these KUCs, we have
also studied the uniaxial compression of a 3 × 2 × 1 kagomé lattice (KL) made from the HKUC
of dimension 70𝑛𝑚 × 70𝑛𝑚 × 70 𝑛𝑚 where each strut has the outer diameter 20 nm and the
thickness 5 nm. The overall size of this KL is 2100𝑛𝑚 × 1400𝑛𝑚 × 70 𝑛𝑚 , and it consists of
30 million atoms. The simulation schedule for the compressive loading of the KL follows the same
schedule as before, and it is compressed up to 15% during the loading phase. Figure 3.17 shows a
schematic of the KL system.
54

Figure 3.17: Schematic of 3 × 2 × 1 kagomé lattice (KL).
3.4.2 Stress-strain curve
Figure 3.18 shows the stress-strain curve of the uniaxial compression of KUCs made from
either solid Ni nanorod or hollow tubes. We will refer the KUC made from the solid Ni nanorods
as SKUC, and the system made from the hollow tubes as HKUC. The Young’s modules of the
SKUC and HKUC is 54.8 GPa and 34.57 GPa, respectively, which is computed from the elastic
regime of the stress-strain curve. The SKUC have higher Young’s modules because it has higher
density and thus have higher load bearing capacity. However, from the stress-strain curve, we can
also observe that yield point of HKUC is 3.9% strain that is higher than that of SKUC’s yield point
of 2.5% strain. Also, both of these systems show differences not only in their Young’s modules
and yield point but also in their response to compressive load after yielding. In case of the SKUC,
load bearing capacity of the system decreases with increasing strain after the yield point. However,
in case of HKUC, the maximum load bearing capacity of the system remains constant up to 11 %
strain and then starts decreasing after that. The difference in their stress-strain curve shows that
55
the HKUC have higher strength and toughness. We have further analyzed the plastic deformation
of these system during yielding and at higher strain in the subsequent sections to grain better
understanding of their deformation mechanism.

Figure 3.18: Stress-Strain of flat punch compression of SKUC and HKUC.
3.4.3 Solid and hollow kagomé cell deformation at yield point
We have analyzed the deformation mechanism of the SKUC and the HKUC by monitoring
the atomic stress and the plastic deformation of the system during the compressive and relaxation
phase of the yield point. Here, plastic deformation is characterized by the nucleation of dislocations
and their density inside the system. To identify the dislocations, we have used the dislocation
extraction algorithm developed by Stukowaski at. el, which characterizes different types of
dislocation present in FCC crystal and identify their burgers vector by computing the centro-
symmetry parameter of each atom.
56

Figure 3.19: Plastic deformation of SKUC at yield point.
In case of SKUC, the deformation of system happens due to the nucleation of
1
6
[112]
Shockley partial and
1
2
[110] dislocation loops on <111> plane inside the system from the
interfacial regions of the KUC with the flat punch and the substrate. These nucleated partial
dislocation loops propagate from the top and bottom interface of the KUC towards its nodal region
through the beams of the kagomé cell on <111> plane (see Figure 3.19). Here, Figure 3.19a
shows the nucleation of a dislocation near the substrate, and Figure 3.19b shows the timestep when
it reaches the nodal region of the KUC. Also, we observe from Figure 3.19b that the new
dislocation loops are nucleated at this time near the flat punch-SKUC interface, which again travels
towards the nodal area of the kagomé cell. At the nodal region, these dislocation loops interact
with each other and some these dislocation loops leaves the system through the surface, which
results into the tiny drop of dislocation density inside the system. However, most of these
dislocation remains inside the system near the nodal region, which now consists of mobile and
immobile partial dislocation loops (Figure 3.19c). These immobile dislocation loops form due to
the interaction between mobile partial dislocation loops. Thus, we observe that nodal area of
kagomé cell is a region of high dislocation density, and hence they also act as a source for
dislocation nucleation along the interfaces. Since, dislocation density does not drop down to zero
(Figure 3.20a), the system does not require higher stress for further plastic deformation at higher
57
strain as observed from the stress-strain curve. Also, we observe that entire SKUC system - the
struts and the nodal region - deforms during the onset of plastic deformation in solid kagomé cell.

Figure 3.20: Dislocation density inside the (a) SKUC and (b) HKUC during yielding.
In case of hollow kagomé cell, the deformation is localized near the nodal area of the
kagomé cell. Here, partial dislocation loops nucleate from the corners of the nodal area and travels
towards the interfaces of the kagomé cell with the flat punch and the subtract. The propagation of
partial dislocation loops follows the similar spiral motion on <111> plane that is also observed
during the mechanical collapse of Ni nanotube using flat punch compression. This happens due to
geometrical constraint of HKUC as it does not have a direct path on the <111> plane from the
nodal region to the interfaces. Thus, the dislocations have to follow a spiral motion, which implies
that they have to travel a larger distance to reach the interfaces. Also, because of this geometrical
constraint, dislocations do not really reach the interface and the deformation is localized only near
the nodal area of the HKUC.
Figure 3.20b shows the dislocation density in the HKUC during yielding where the drop-
in dislocation density happens due to the annihilation of dislocation loops from the inner and outer
58
walls of the kagomé cell, however, the dislocation density does not drop down to zero just like the
SKUC and here also higher dislocation density is observed in the nodal region. In contrast to solid
kagomé cell, only the nodal region of the hollow kagomé cell deforms, and there is not much
deformation in of its struts.
From the stress analysis in these systems before yielding, we observe that solid kagomé
cell has higher critical resolved shear stress near the nodal region and its interfacial region with
the flat punch and substrate; however, in case of hollow kagomé cell maximum stress is near its
nodal area. Figure 3.21 shows regions of the SKUC and HKUC with stress greater than 3GPa.

Figure 3.21: Shows the regions of stress distribution with stress greater then 3GPa inside the (a)
solid kagomé cell and (b) hollow kagomé cell just before yielding.
3.4.4 Local deformation analysis
The analysis of plastic deformation of the systems by dislocation activity requires us to
analyze all the frames so to get an estimate of the locations of deformed regions inside the system.
However, we still cannot visualize all the deformed regions in a single frame. To find these
locations of plastic deformation by just analyzing a single frame, we can use local deformation
analysis which is developed by M.L. Falk. Local deformation analysis can provide quantitative
measure of the locations inside the system, which is affected by the plastic deformation by just
looking at the last frame at the end of a compressive loading. It computes local deformation
59
parameter D
2
min for each atom by using the initial and final coordinates of atoms inside the system
before and after the compressive loading, which can then be used to quantify deformed regions
inside the material. D
2
min for each atom is calculated by first defining an interaction range (rc) for
each atom, which here we have taken 6 Å. After that, local strain, 𝜖 𝑖𝑗
, around each atom is
calculated by minimizing the mean square difference between the displacements of the center atom
(i) and its neighbors inside the interaction range in the deformed system and the relative
displacements that they would have if they were under uniform strain 𝜖 𝑖𝑗
in some other reference
frame. Here the initial system before the application of compressive load is taken as the reference
frame. D
2
min for each atom is calculated using the equation 1
𝐷 2
(𝑡 ,∆𝑡 )=∑ ∑ {𝑟 𝑛 𝑖 (𝑡 )−𝑟 0
𝑖 (𝑡 )−∑ [(𝛿 𝑖𝑗
+𝜖 𝑖𝑗
)∗(𝑟 𝑛 𝑖 (𝑡 −∆𝑡 )−𝑟 0
𝑖 (𝑡 −∆𝑡 ))]
3
𝑗 =1
}
2
3
𝑖 =1
𝑁 𝑛 =1
Eq. 8
Here, r0 is the coordinate of the atoms in the reference frame and rn is their coordinate in
the deformed system at time t and the minimum value of 𝜖 𝑖𝑗
is determined by equation 2-4.
𝑋 𝑖𝑗
=∑ [𝑟 𝑛 𝑖 (𝑡 )−𝑟 0
𝑖 (𝑡 )]
𝑛 ∗[𝑟 𝑛 𝑗 (𝑡 −∆𝑡 )−𝑟 0
𝑗 (𝑡 −∆𝑡 )] Eq. 9
𝑌 𝑖𝑗
=∑ [𝑟 𝑛 𝑗 (𝑡 )−𝑟 0
𝑗 (𝑡 )]
𝑛 ∗[𝑟 𝑛 𝑗 (𝑡 −∆𝑡 )−𝑟 0
𝑗 (𝑡 −∆𝑡 )] Eq. 10
𝜖 𝑖𝑗
=∑ 𝑋 𝑖𝑘
𝑌 𝑖𝑘
−1
−𝛿 𝑖𝑗 𝑘 Eq. 11
D
2
min gives us the local deviation from affine deformation and is used as measure of plastic
deformation in the system. Figure 3.22 shows the regions of the system in solid and hollow kagomé
cell where D
2
min is greater than 0.25 after yielding. We can observe that local deformation analysis
is able to locate the regions of high plastic deformation inside the system by just analyzing the last
frame after yielding. These regions are highly localized near the nodal region in hollow kagomé
cell, and in case of solid kagomé cell the deformed regions are in the struts and nodal area. It also
60
tells us the different planes in the system where the deformation has happened, which are here
system of {111} planes and their intersections.

Figure 3.22: Local deformation analysis of (a) SKUC and (b) HKUC showing regions of plastic
deformation after yielding.
Local deformation analysis is very useful for large scale simulations involving millions to
billions of atoms as it does not require us to save all the frames during the simulation for post
processing and data analysis.
3.4.5 Solid and hollow kagomé cell deformation at higher strain
Figure 3.23 shows local deformation analysis of solid kagomé cell for compressive strain
at 12% and 19% strain. We can observe from figure 3.23 that large regions inside the system
suffers from plastic deformation at higher strain. Here, deformation happens again inside all the
eight struts and the nodal area of the solid kagomé cell. Due to the large deformation inside the
system, we observe formation of several types of defects like twinned and slipped regions inside
the material as shown in figure 3.23b and 3.23d.

61

Figure 3.23: Local deformation analysis showing regions of plastic deformation in SKUC at 12%
and 19% strain. (a) and (c) shows internal regions with D
2
min greater than 0.25, and (b) and (d)
shows outer surface of SKUC, which shows twin boundaries.
Figure 3.24 shows local deformation analysis of hollow kagomé cell for compressive strain
at 11% and 19% strain. Unlike SKUC, KHUC’s deformation follows a two-stage deformation
mechanism where up to 11% applied strain, deformation is localized near the nodal region. From
the stress-strain curve we can observe that up to 11% strain, maximum load bearing capacity of
the hollow kagomé cell remains almost constant, and this happens due to this highly localized and
little deformation near the nodal area. At strain above 11%, we observe both the nodal region and
the eight struts start deforming. We observe large deformed regions inside the system. Because of
these large plastic deformations inside the system, we observe drop in load in the stress-strain
curve for strain above 11%. Stress analysis in the hollow kagomé cell also shows that up to 11%
strain, maximum stress inside the system is near the nodal region while at strain above 11% it’s
near the node and the interfaces.
62

Figure 3.24: Local deformation analysis showing regions of plastic deformation in HKUC as a
function of applied strain.
3.4.6 Comparison of kagomé cell deformation with nanorod deformation
Figure 3.25 shows the dislocation density inside the solid and hollow kagomé cell as a
function of applied strain. We can observe that dislocation density inside these systems increases
with the increase in applied strain, which is very different from the deformation mechanism
observed in nanorods and nono-tubes. In nanorod and nano-tube, dislocation nucleation inside the
system follows dislocation starvation mechanism. Dislocation starvation mechanism says that due
to the smaller length of the nanorod and nanotube, during deformation, dislocation nucleates from
the surfaces of the nanotube/nanorod and travels through the material. However, before these
dislocations does not have enough time to interact with other dislocations and multiply, they leave
the material again from the surface. Hence, the system remains in a dislocation starved state after
the deformation, and it again higher stress for the nucleation of new dislocations. In solid and
63
hollow kagomé cell, from the dislocation density plot (Figure 3.25), we can observe that
dislocation density inside the material never drops down to zero at any strain, in fact it increases
with applied strain. Also, in both system, maximum dislocation density is localized near the nodal
regions. These dislocations can interact and create new dislocation at higher strain along with the
dislocation nucleation from the surface of kagomé cell.

Figure 3.25: Dislocation density inside the SKUC and HKUC as a function of applied strain.
3.4.7 Bending of solid and hollow kagomé cell
Dominant deformation mechanism observed in hierarchical material are either stretching-
dominated or bending-dominated deformation. Here, we have observed bending of the struts of
kagomé cell in response to the compressive load. Similar bending dominated deformation is
observed during the flat punch compression of kagomé lattice made from ploy-crystalline Ni by
Schaedler et. al. Hence, we have further characterized the deformation of the solid and hollow
kagomé cell by the amount bending of its struts in response to compressive load. Figure 3.26 shows
the bending of the beams of the kagomé cell as a function of applied strain. We observe that solid
kagomé cell shows higher bending. In case of hollow kagomé cell, the beams of the kagomé cell
64
shows very little bending up to 11% strain and it increases rapidly after that. However, the overall
bending observed in the system is still smaller than that of solid kagomé cell. Smaller bending of
the hollow kagomé cell’s beams up to 11% strain is due to the highly localized deformation near
the nodal region and very little deformation in its beams. These smaller bending and plastic
deformation of the hollow kagomé cell shows that it has higher strength and toughness in
comparison of solid kagomé cell.

Figure 3.26: Angle 𝜃 between the struts at the node as a function of applied strain in SKUC and
HKUC.
3.4.8 Uniaxial compression of hollow kagomé lattice
Figure 3.27(a) shows the stress-strain of the uniaxial compression of the hollow kagomé
lattice, figure 3.27(b) shows the system after 15% deformation. From stress-strain curve, we
observe that after initial yielding, the stress first decreases and then settles down at a constant value
of 0.8 GPa. Similar stress-strain curve is observed experimentally for other bending dominated
micro lattices. This drop happens due to the onset of plastic bucking in the system.
65

Figure 3.27: (a) Stress-Strain curve of kagomé lattice, (b) regions of plastic deformation inside
the KL at 15% strain. Here, green areas are stacking fault and red lines are dislocations.
From Figure 3.27, we can observe that, most of the deformation in the system is localized
near the middle kagomé cell while very little deformation has happened in top and bottom cells.
3.4.9 Summary
In summary, we observe that hollow kagomé unit cell (3.9 % strain) has higher yield point
then solid kagomé unit cell (2.5 % strain). HKUCs shows two stage deformation process where
up to 11% strain, the deformation is highly localized near the nodal region and for strain above
11%, all the struts and node starts deforming. Due to this two-stage deformation mechanism
HKUCs shows less bending and higher toughness in comparison of SKUCs.
3.5 Conclusion
Molecular dynamics simulation shows that deformation mechanism of nanorod and
nanotube is very different from each other. As compared to nanorod, in nanotube deformation
happens due to spatial motion of partial dislocation loops on (111) plane due to its geometrical
66
constraints. Because of this unique dislocation activity in nanotube, the size of the deformation
zone (twin regions, stacking fault) in nanotubes are smaller in comparison of nanorods. The
nanorod shows mechanical collapses due to the nucleation of partial dislocations loops at from the
nanorod surface at 3.5% strain.
This difference in deformation mechanism of the nanorod and the nanotubes makes the
response of SKUC and HKUC contrastingly different. HKUCs shows a two-stage deformation
process where up to 11% strain, the deformation is highly localized near the nodal region and for
strain above 11%, all the struts and node starts deforming. Due to this two-stage deformation
mechanism HKUCs shows less bending and higher toughness in comparison of SKUCs.
A KL lattice made from HKUC shows even unique deformation mechanism where system
is able to withstand up to 15% with localized deformation in the middle kagome cells and very
little deformations in the entire system.

67
Chapter 4: Self-Healing Ceramics
4.1 Background
There is a great deal of interest in designing material systems having the ability to sense and
repair damage.
1,10,52-54
Usually damage initiation involves the formation of defects such as
dislocations, voids and micro-cracks. These defects grow and coalesce to cause material failure. If
we can manually detect and repair these micro-cracks at early stage, then it will enhance their
service life, which in many situations is not a feasible solution because it requires constant human
intervention. In this respect, we can learn a lot from nature as many natural materials shows
automatic healing phenomena for example, plants and animals have the ability to heal their
wounds.
9
It is not unreasonable to say that many materials (metals, ceramics) have superior
mechanical properties (toughness, strength) in comparison of the biological materials like bone,
tendon, arthropod cuticle, plant and stem shells. Yet natural materials have high structural
integrity, can withstand high stress and have longer service life. Biological materials have two
unique characteristics that help them to achieve these properties. These are hierarchy in their
micro-structure and self-healing phenomena.
There are many examples of micro-structural hierarchy in biological material that help
them achieve high strength and toughness.
9
For example, shell of sea creatures (e.g. mussel, limpet
and conch, nacre) is made from calcium carbonate (ceramic) with small amount of water and
organic materials. Even though calcium carbonate has very low fracture toughness (KIC= 0.3),
which is similar to that of egg shells, but nacre has higher toughness. This happens due to the
complex hierarchy of their microstructure that consists of tile-like units at the micron scale with
specially developed interface. Another example is human bone whose multi-level hierarchical
structure enhances their toughness via crack defection and crack-bridging ligaments. This micro-
68
structural hierarchy in material have been utilized by materials scientists to design ultra-low
density and high strength materials.
The second feature observed in several biological materials is their self-healing ability that
further increases their reliability and damage tolerance. For example, the load bearing tissues in
our body like bones and skins shows continuous healing at microscopic level which helps tour
body to reduce long-term failure mode like fatigue and fracture. Taking this cue from biology,
material scientists are developing novel self-healing systems to enhance reliability and lifetime of
materials while reducing the cost of manufacturing, monitoring and maintenance.
1,52-54
Most of the
approaches to designing self-healing material systems involve nanoscale or mesoscale containers
filled with a healing material, which is released upon damage initiation either by an external
stimulus (heat, light, electric or magnetic field) or an internal change in the state of the system (pH,
stress or temperature).
11-14
The healing material is transported to the damage site and immobilized
after healing damage. Different materials require different self-healing conditions — ambient
conditions for concrete, low temperatures for polymers, relatively high temperatures for metals,
and very high temperatures for ceramics.
Self-healing polymers
Polymers and polymer matrix composites (PMCs) have widespread application in many
areas like aerospace, consumer goods, packaging and as structural material
10,11,55
. One such
example is the Boeing 787 airplane, which is completely made from the PMCs. In general,
common failure modes observed in these PMCs are delamination, fiber-matrix debonding, fiber
fracture and micro-cracking of the brittle polymer matrix, which limits their full potential to be
realized. Conventional approach for failure prevention in these materials includes damage
detection, damage prediction, protective coating and manual repair but all these techniques are
69
time consuming, expensive and require human intervention. PMCs that mimic biological materials
is another design choice as they can automatically heal minor cracks and voids, therefore increase
their service life. Since early 1990’s, self-healing mechanism in polymers and PMCs have been
extensively studied by material scientists and can be categorized into four groups; (1) healing by
crack-filling adhesion, (2) healing by diffusion, (3) healing by bond reformation and (4) virgin
property strengthening in response to stress.
Crack-filling adhesion
In crack-filling adhesion, a liquid substance is released upon damage initiation, which
travels to the damage region and completely fills. Here, crack healing happens due to the various
chemical and physical processes which occurs in these liquid substances once they come in contact
with the crack surface. The healing substance used here includes liquid monomers, liquid catalyst,
thermoplastic polymers, organic and inorganic film formers that are collectively called healing
agent. In comparison of other healing mechanism observed in polymers, it does not require crack
surfaces to come in contact for healing to initiate as here repair happens due to the flow of healing
agent to the damaged region. It can be further divided into four types, which are microencapsulated
healing agent based, phase separated healing agent based, healing agent filled inside hollow fibers
and microvascular network.
In microencapsulated based healing (Figure 4.1a), microcapsules containing liquid healing
agent and catalysts (either solid particles or another encapsulated liquid) are randomly distributed
throughout the polymer matrix. When cracks start propagating through the material, they rupture
many microcapsules that come across their path. Due to capillary actions, the catalyst and the
healing agents inside these ruptured microcapsules goes the crack where the healing agent
polymerizes in the presence of catalyst and adheres the two crack faces together. Here, it is
70
important that the healing agent and microcapsules remain stable during the processing of the
composite. Further, the healing agent should have low viscosity so that it can travel faster to the
damage sites upon its release. An example of such material is first demonstrated by White in 2001,
where crack healing in epoxy matrix is shown using a liquid healing agent dicyclopentadiene
(DCPD) encapsulated inside the poly(urea-formaldehyde) shell. They have also used Grubbs’
catalyst, which is dispersed separately along with the microcapsules throughout the epoxy matrix.
Another example demonstrated by Rong et al. was that both the healing agent and the base material
were the same polymer so as to better mimic the biological materials. They have used a diglycidly
ether bisphenol-A (DGEBA) based epoxy resin encapsulated inside urea-formaldehyde
microcapsule along with imidazole hardener embedded inside epoxy resin, which is made from
the same DGEBA. These systems showed 100% recovery of fracture toughness after healing. The
healing temperature required for these materials after damage is between room temperature to 140
°C where higher curing temperature is required when imidazole hardener is used. The curing time
for crack healing is between 8 to 15 hours, which in general depends upon the healing agent.
In phase separated healing agent-based polymers, healing agent is directly incorporated
inside the polymer matrix as a separate phase. In comparison of microcapsuled based healing, they
have relatively simpler processing step and requires less fabrication time as the encapsulation step
is not necessary. However, the healing agent has to be inert, stable and should not reach with the
matrix, which is important as they do not have protective encapsulation. An example of such
material is first demonstrated by Zako et al. They have created a fiber-reinforced polymer matrix
containing nanoparticles of the epoxy prepolymer as healing agent. To mitigate the damage regions
in the matrix, it is heated above the melting temperature of the healing agent (~120 °C), where
71
epoxy prepolymer on the fracture surface melts, fills the damaged regions. Upon cooling these
epoxy prepolymer cures and heals the polymer matrix with complete recovery of strength.

Figure 4.1: Self-healing phenomena in polymers (a) microencapsulated based healing, (b) hollow-
pipes containing healing agent and (c) microvascular network-based crack healing
55

The above-mentioned techniques do not allow repeated healing of a specific area in the
polymer matrix because healing agent inside the microcapsules have limited quantity, which once
depleted cannot provide any more recovery to the PMCs. To handle this situation, material
scientists have designed polymer composite embedded with hollow fibers. These hollow fibers act
as reinforcement and also contain liquid healing agent. Upon damage, they are able to supply
healing agent not only to a large localized damaged region but also continuously to a specific
region for a longer period of time. Healing agent is introduced into the hollow fibers through its
open ends, capillary action, vacuum assistance or through surface pores, which is covered later. A
schematic of such composite is shown in Figure 4.1b. A polymer composite containing hollow
borosilicate glass fibers is developed by Bond et at. These glass fibers were filled with epoxy
resins and UV fluorescent dye as healing agent, and they are aligned as groups of orthogonal pipes
inside the matrix so as to provide strength.
Based on the concept of vascular system of plants and animals, material scientist has
designed polymer composite that contains interconnected network of hollow pipes, which is
72
known as microvascular network. These microvascular network mimics vascular network of
humans and are able to provide continues supply of healing agent to all the location inside the
matrix. To realize such micro supply-network, they can be connected to an external, refillable
liquid pump that can deliver a constant supply of healing agent. Hence, they are superior in
comparison of all other technique in which healing agent eventually get exhausted. Figure 4.1c
shows a schematic of such network in a polymer matrix. However, design of such microvascular
network architecture is a complex task as it needs to be optimized so that healing agent can be
delivered to multiple damage site with maximum flow properties. A fabrication technique of such
material is based on first developing a scaffold of fugitive wax, which is then infused with an
epoxy/hardener resin. Later the wax is removed by heating it to higher temperature, which creates
microchannels inside the epoxy matrix. These microchannel are then filled with healing agent
(DCPD). Another design technique consists of creating a PMCs where PVC tubes are embedded
horizontally inside the matrix. Here, microvascular system is created by drilling vertical holes
inside these tubes at specific locations.
Crack healing by diffusion
In these systems, crack healing happens due to the diffusion of healing agent to the damage
where they form chemical and physical linkage to do crack closure. Unlike crack filling adhesion
where healing happens due to the flow of healing agent to the damage site, this process requires
external stimulus for the transport of healing argent from one damage surface to another. For
example, polymer matrix is created such that it has enough dangling bond which can be achieved
by eliminating the sol phase of weak polyurethane sol-gel just below its critical point. When
cracked surface of such polymers are bought into contact then due to the topological interaction
between dangling chains, system shows crack closure. Another example of such material is shown
73
by Rahmathullah and Palmese using thermos epoxy amine polymers. In these materials, when
damaged surfaces are brought into contact under low pressure and at temperature above 185
0
C,
they show crack closure due to localized softening. These localized softening increases molecular
mobility at the crack surface and hence helps in crack closure. Also, some materials utilize their
viscoelastic property for crack closure. When these materials are subjected to high energy impact,
the energy transferred to the material due to the impact causes localized melting and helps in crack
closure.
Crack healing by bond reformation
Some polymers show unique characteristics as they contain special bonds that are
reversible under extrenal stimulus. These polymers are called monomers. They are utilized to
design materials in which under stress the reversible bonds of monomers break, which is recovered
again using external stimulus, and the material regains its strength. Various external stimulus used
for crack healing are thermal stimulus, electrical signal and ultra-violet light. For examples,
healing property is shown by materials made from diene and dienophile who reacts via Diel-Alder
(DA) reaction. However, their reaction product acts as a weak link during damage, and the DA
adduct reverts back to the original molecules diene and dienophile. Here, healing happens by DA
reaction again in presence of heat at around 120 °C. Murphy et al. have sensitized self-healing
monomers based on this reaction, where derivative of cyclopentadiene is used as diene and
dienophile. Examples of materials in which UV-light is used as an external stimulus is polymer
materials containing tricinnamters. These tricinnamters molecules form cyclobutene, which
contains cross-linked polymer films. These cross-linked polymer films break down to original
tricinnamters during crack propagation, but the crosslinking is restored again in the presence of
UV-light.
74
Virgin property strengthen
In these class of polyremes, strengthening happens due to the localized polymerization
reaction and cross-linking in the damaged region in response to applied stress but before virgin
material failure. These polymers contain mechanophores, which are mechanically active units and
are presents along the backbone of the polymer. Examples include ruthenium bis-carbene
complexes with polymer functionalized legends which shows strengthening via stress induced
olefin metathesis reaction.
Self-healing phenomena in metals
In comparison of polymers, self-healing phenomena in metals and ceramics is often more
complex and has not been extensively studied because of their high operating temperature and high
melting temperature resulting in the low atomic mobility. Since metals operate at high temperature,
they generally fail via low-ductility creep fracture. These happens due to the nucleation and growth
of voids that coalesce and form cracks and results into creep failure. To prevent high-temperature
creep damage in metals, material scientists use precipitate hardening of the matrix, where healing
happens due to the precipitation of solute particles on the free surfaces of voids. These precipitates
fill the void completely and stops crack growth. For example, underaged aluminum alloy shows
higher creep resistance in comparison of fully hardened aluminum alloy at 150 °C.
56
This happens
due to enhanced dynamic precipitation of excess residual solute, which reduces secondary creep
deformation. Another example of higher creep resistance is observed in iron containing small
amount of gold.
57
These mobile gold particles inside the iron matrix increases its service life at
temperatures of up to 550 °C. This happens due to selective precipitation of gold particles on free
surfaces of cavities, and hence acts as autonomous repair of the creep damage. Here grain
75
boundaries, dislocations act as fast medium for the transportation of the gold particles to the
damage site.
Self-healing phenomena in ceramics
Ceramics are lightweight, have high strength and stiffness, and can withstand high
temperatures. However, technological applications of ceramics are limited by brittleness and
sensitivity to flaws and cracks. Self-healing of defects and cracks can dramatically increase the
reliability and lifetime of ceramics which, in turn, can reduce maintenance costs for a broad range
of energy technologies from ceramic turbine blades to solid-oxide fuel cells.
58,59
Strength recovery
in ceramics is categorized into three groups that are: (1) re-sintering, (2) relaxation of residual
stress and (3) crack bonding by oxidation. Re-sintering of ceramic components decreases the flaw
size but sometimes it generated spherical voids inside the material. Relaxation of residual stress is
usually carried out in vacuum, air, nitrogen and argon atmosphere. It increases the strength of the
material, but it does not heal the cracks insides the material and thus it can still suffer from
catastrophic failure. Crack bonding by oxidation has been studied extensively in various ceramic,
which shows that they recover their strength after this treatment. Examples of crack healing by
oxidation in ceramics includes mullite-SiC composite, Si3N4-SiC composite,
16,17
Al2O3-SiC
composite,
18,19
MAX phase ceramics (Ti3AlC2, Ti2AlC)
20,21
and also in an aluminosilicate glass
containing vanadium-boride nanoparticles.
60,61

All ceramics that shows oxidation-based crack healing, contains two ceramic phases where
one phase reacts with the oxygen from the air and forms oxide. These oxides fill the crack and
heals the system completely. For example, in mullite-SiC, Si3N4-SiC, Al2O3-SiC, SiC reacts with
the oxygen at higher temperature and forms molten silica who fills the crack completely. In
general, SiC is dispersed as nanoparticle or whiskers inside the ceramic matrix. Figure 4.2 shows
76
schematic of SiC oxidation-based crack healing and the governing equation of this reaction is
shown in equation 4.1.
𝑆𝑖𝐶 (𝑠 )+
3
2
𝑂 2
(𝑔 )→𝑆𝑖 𝑂 2
(𝑠 )+𝐶𝑂 (𝑔 ) Eq4. 1

Figure 4.2: Crack healing due to flow of liquid silica into the surface cracks.
In MAX phase ceramics, crack healing happens due to the formation of Al2O3 by oxidation
in air at 1200 °C. For example, in Ti2AlC, oxidation happens due to inward diffusion of oxygen
inside the matrix and outward diffusion of Al and Ti by following reaction:
4𝑇 𝑖 2
𝐴𝑙𝐶 (𝑠 )+15𝑂 2
(𝑔 )→8𝑇𝑖 𝑂 2
+2𝐴 𝑙 2
𝑂 3
+4 𝐶 𝑜 2
Eq4. 2
Ceramic materials develop surfaces crack during their manufacturing, which severely
reduces their service life. To increase the reliability of these ceramic composite three step
procedure is applied. In step 1, surface cracks are automatically healed using crack healing
treatment at higher treatment (~1,200 C). In second step using proof-test, components with non-
acceptable flaw are discarded. In the last step, if crack develops during their service life of the
material, they should be able to automatically heal those cracks too. To do that, experiments has
been done to determine the minimum threshold stress that the material can withstand and show
crack healing. For example, Ando et al. have shown the effect of stress on crack healing in
77
Al2O3/SiC composite containing indentation crack.
62-65
For this, they have done three-point
bending measurement of the samples at 1,300 C as a function of stress applied on the samples.
They have shown that if the applied stress is less than a critical value of 250MPa, it shows crack
healing phenomena. This threshold stress depends on the material and is given by Irwin criteria.
According to Irwin, the driving force for crack growth increases with crack length and applied
stress. But the driving force for crack healing is independent of applied stress, and it increases with
healing temperature. Hence crack healing is able to prevent crack growth below a critical applied
stress at temperature regime where crack healing rate dominates the crack growth rate.
In general crack healing depends on five parameters, which are (1) composition of the ceramics,
(2) SiC particle size, (3) atmosphere, (4) temperature and (5) applied stress on the system.
15
It is
observed that SiC shape, size and its composition are very important for crack healing. For
example, Ando et al. found that the critical concentration of n-SiC for crack healing in Al2O3 is
between 15- 20%. Experiments by Niihara et al. indicate that five percent volume fraction of sub-
micron size SiC particles can increase the fracture toughness of alumina by 40% and strength by
1 GPa.
66
Similarly spherical particles are more effective in crack healing in comparison of
whiskers. This happens because whiskers forms at the crack wall and reduces the effective
oxidation at the crack surface. It is observed that crack healing happens maximum in air as it
contains oxygen. Other environments like N2, vacuum, Ar are not very effective. Also, the required
temperature for crack healing follows the Arrhenius equation, and required time for crack healing
decreases with temperature.
Crack healing in Al2O3
Al2O3 has excellent mechanical properties (hardness, wear resistance) and oxidation resistance
compared to metals at elevated temperature. However, it has following weakness low strength
78
(  = 400Mpa), low fracture toughness (KIC = 3MPa m
1/2
) and low heat resistance limit temperature
for strength (1200K). Ando et al. has shown experimentally that these weaknesses of Al2O3 can
be alleviated by creating nanocomposite of Al2O3 containing SiC nanoparticles.
62-65
These
Al2O3/SiC nanocomposite shows crack healing property and because of that they have superior
mechanical property. In their experiment, Ando et al. introduced indentation cracks on Al2O3
surface and preformed three-point measurements to demonstrate the recovery of strength in crack-
healed samples. After creating the indentation crack, the samples are heated up to 1,300 C for 2
hours so as to heal the indentation crack. On these crack healed samples, they have performed
three-point bending experiment as a function of temperature from 25 C to 1,500 C. Their
experiments have shown that bending strength of Al2O3 decreases with increasing temperature,
and the deformation in Al2O3 is plastic in nature above 1,400 C. Also, in all these crack healed
samples, fracture initiated from outside the crack-healed zone compared to the systems which did
not go through these crack healing treatment and fracture happened in these samples just below
the indentation crack. Figure 4.3 shows the fracture of a composite after crack healing treatment.
They have also shown that the concentration of SiC nanoparticle inside Al2O3 is critical for crack
healing and it should be at least 15-20% by volume fraction. Experiments by Niihara et al. indicate
that five percent volume fraction of sub-micron size SiC particles can increase the fracture
toughness of alumina by 40% and strength by 1 GPa.
66

79

Figure 4.3: (a) Indentation crack on Al2O3 surface before crack healing treatment. (b) Crack
healed sample showing fracture at the different location than the original indentation crack.
The suggested mechanism for crack-healing here is the formation of silica shells around n-
SiC by oxidation of Si with atmospheric oxygen followed by diffusion of liquid silica into cracks.
The exothermic heat and volume expansion associated with oxidation of n-SiC are believed to play
important roles in accelerating crack healing by diffusion of silica into the damaged zone.
However atomistic mechanism of crack healing in Al2O3 is not completely clear. Hence,
here we examined the atomistic mechanisms of crack healing in Al2O3 by oxidized n-SiC and the
effect of nanoparticle size and distribution on crack healing.

80
4.2 Simulation details
4.2.1 Crystal structure of Al2O3 and SiC
-Al2O3 has hexagonal crystal structure with the unit cell dimensions a = 0.4754 nm and
c = 1.299 nm. Transforming the hexagonal unit cell into an orthorhombic cell, we obtain the
orthogonal unit cell dimensions of -Al2O3 will be a = 0.4754 nm, 𝑏 =√3𝑎 =0.8234 𝑛 𝑚 and c
= 1.299 nm. In the -Al2O3 crystal, the atomic coordinate of the anions forms hexagonal sub-
lattice along the basal plane and hence have an …ABABAB… stacking whereas cations take 2/3
of the octahedral interstitial sites along the basal plane and form an …ABCABCABC… stacking.
Figure 4.4a shows this stacking of anions and cations along the basal plane, where the viewing
direction is the prism plane of Al2O3. Figure 4.4b shows the atomic arrangement when the viewing
direction is along the basal plane. In this direction, atoms show hexagonal arrangement.

Figure 4.4: −Al2O3 crystal along the (a) prism plane and (b) basal plane.
SiC has many different crystal polymorphs, for example 3C-SiC (Zinc blende, cubic), 4H-
SiC (Wurtzite, Hexagonal), 6H-SiC (Wurtzite, Hexagonal), and 15R-SiC (Rhombohedral). Here,
81
we have used 3C-SiC crystal structure whose space group is 𝑇 𝑑 2
−𝐹 4
̅
3𝑚 and lattice constant is
0.43596 nm.
4.2.2 Interaction potential for Al2O3, SiC and SiO2
We have used the reactive force field developed by Vashishta et al. to simulate the crack
healing mechanism in an -Al2O3 matrix containing silica-coated SiC nanoparticles.
67-69
This
reactive force field consists of a two-body and a three-body term to incorporate various effects,
including steric repulsion Coulomb interaction, charge-dipole interaction, van der Waals
interaction, covalent bond bending and stretching, and it allows bond breaking and bond formation.
The functional form of this interaction potential is written as:
𝐸 𝑡𝑜𝑡 =∑ 𝑉 𝑖𝑗
(2)
(𝑟 𝑖𝑗
)+∑ 𝑉 𝑖𝑗𝑘 (3)
𝑖 <𝑗 <𝑘 𝑖 <𝑗 (𝑟 𝑖𝑗
⃗⃗⃗ 𝑟 𝑖𝑘
⃗⃗⃗⃗ ) Eq4. 3
Where, 𝐸 𝑡𝑜𝑡 is the total potential energy of the system, and 𝑟 𝑖𝑗
⃗⃗⃗ =𝑟 𝑖 ⃗⃗ −𝑟 𝑗⃗⃗ is the separation
between atoms i and j positioned at 𝑟 ̅
𝑖 and 𝑟 ̅
𝑗 , respectively. The first term in equation 4.3 is the two-
body terms, and it consists of steric repulsion, screened Coulomb, charge-dipole, and van der
Waals 1/r
6
terms. The potential energy of an atom, i, due to this 2-body term, sums over all its
neighbors j that are within the cutoff distance rc from its position. Its functional form is given in
equation 4.4.
𝑉 𝑖𝑗
2
(𝑟 )=
𝐻 𝑖𝑗 𝑟 𝜂 𝑖𝑗 +
𝑍 𝑖 𝑍 𝑗 𝑟 𝑒 −
𝑟 𝑟 𝑖𝑠 −
𝐷 𝑖𝑗 2𝑟 4
𝑒 −
𝑟 𝑟 4𝑠 −
𝑤 𝑖𝑗 𝑟 6
Eq4. 4
The second term in equation 4.3 is the contributions due to the three-body terms. Here, in
the three-body interaction, we sum over all the atoms, j and k, that are covalently bonded with
atom i within a distance r0, and is written as:
𝑉 𝑖𝑗𝑘 (3)
(𝑟 𝑖𝑗
⃗⃗⃗ 𝑟 𝑖𝑘
⃗⃗⃗⃗ )=𝐵 𝑖𝑗𝑘 exp(
𝜉 𝑟 𝑖𝑗
−𝑟 0
+
𝜉 𝑟 𝑖𝑘
−𝑟 0
)
(𝑐𝑜𝑠 𝜃 𝑖𝑗𝑘 −𝑐𝑜𝑠 𝜃 0
)
2
1+𝐶 𝑗𝑖𝑘 (𝑐𝑜𝑠 𝜃 𝑖𝑗𝑘 −𝑐𝑜𝑠 𝜃 0
)
2
(𝑟 𝑖𝑗
,𝑟 𝑖𝑘
≤𝑟 0
) Eq4. 5
82
In equation 4.5, Hij is the steric repulsion strength, Zi is the effective charge in units of the
electronic charge, Dij represents the charge-dipole strength, Wij is the strength of van der Waals
attraction. 𝜂 𝑖𝑗
represents the exponent of the steric repulsion, and r1s and r4s are the screening length
for the Coulomb and charge-dipole interactions, respectively. In equation 4.5, Bijk is the strength
of the three-body interaction, 𝜉 and Cijk are constants. 𝜃 𝑖𝑗𝑘
is the angle between the two vectors 𝑟 𝑖𝑗
⃗⃗⃗
and 𝑟 𝑖𝑘
⃗⃗⃗⃗ spanned by an atomic triplet i, j, k, and 𝜃 0
is their equilibrium angle. For computational
efficiency, the two-body term, 𝑉 𝑖𝑗
2
(𝑟 ) , as shown in equation 4.4 is truncated at the cutoff distance
rc and the contribution from atoms beyond rc is ignored. This truncation makes the derivative of
𝑉 𝑖𝑗
2
(𝑟 ) discontinues at rc. To make the derivative of 𝑉 𝑖𝑗
2
(𝑟 ) at rc, it is also shifted for r < rc and is
written as:
𝑉 𝑖𝑗
2(𝑠 ℎ𝑖𝑓𝑡𝑒𝑑 )
(𝑟 )={
𝑉 𝑖𝑗
2
(𝑟 )−𝑉 𝑖𝑗
2
(𝑟 𝑐 )−(𝑟 −𝑟 𝑐 )
𝑑 𝑉 𝑖𝑗
2
(𝑟 )
𝑑𝑟
|
𝑟 =𝑟 𝑐 (𝑟 <𝑟 𝑐 )
0 (𝑟 ≥𝑟 𝑐 )
Eq4. 6
The parameters of this force field are fitted such that it can correctly describe the various
physical and mechanical properties of SiC, Al2O3 and amorphous SiO2. For Al2O3, the parameters
of the force field were chosen to reproduce the correct lattice constant, cohesive energy, structural
and mechanical properties of crystalline and amorphous phases and melting behavior of -Al2O3.
The parameters of the potential -Al2O3 is given in Table 4.1 and Table 4.3 shows the comparison
of various properties of -Al2O3 with the experimental data.
Al O
Zi (e) 1.5237 -1.0158
r1s = 5.0 Å r4s = 3.75 Å rc = 6.0 Å e = 1.60210
-19
C
Two body
83
Al-Al Al-O O-O
ij 7 9 7
Hij (eV Å

) 12.4236 244.3038 553.3917
Dij (eV Å
4
) 0 3.4374 1.52775
Wij (eV Å
6
) 0 0 63.4894
Three body
Bijk (eV) o (deg) Cijk  (Å) r0 (Å)
O-Al-O 12.4844 90.0 10 1.0 2.90
Al-O-Al 8.1149 109.47 10 1.0 2.90
Table 4.1: Parameters of the Vashishta interaction potential for Al2O3.

Al2O3
Vashishta Potential
Expt. DFT
B (GPa) 253 255 246.9
C11 (GPa) 523 498 476.8
C12 (GPa) 147 164 157.6
C13 (GPa) 129 117 119.4
C14 (GPa) 7.5 -23 19.4
C33 (GPa) 427 502 476.6
C44 (GPa) 135 147 145.5
C66 (GPa) 174 167 159.6
Poison ratio (𝜈 ) 0.22 0.231 0.2356
Table 4.2: Comparison of mechanical properties of Al2O3 predicted by Vashishta potential with
the experimental and DFT data.
84
Si C
Zi (e) 1.201 -1.201
r1s = 5.0 Å r4s = 3.00 Å rc = 7.35 Å e = 1.60210
-19
C
Two body
Si-Si Si-C C-C
ij 7 9 7
Hij (eV Å

) 23.67291 447.09026 471.74538
Dij (eV Å
4
) 2.1636 1.0818 0
Wij (eV Å
6
) 0 61.4694 0
Three body
Bijk (eV) o (deg) Cijk  (Å) r0 (Å)
Si-C-Si 9.003 109.47 5.0 1.0 2.90
C-Si-C 9.003 109.47 5.0 1.0 2.90
Table 4.3: Parameters of the Vashishta interaction potential for SiC.
The parameters of the interatomic potential for SiC were chosen to reproduce the correct
lattice constant, melting point, vibrational density of states, elastic moduli, and cohesive energy of
3C-SiC crystal. Table 4.3 The parameters of the SiC interaction, and Table 4.4 shows the
comparison of various properties of SiC with the experimental data.
3C-SiC Vashishta Potential Expt. DFT
B (GPa) 225.1 225 224.9
C11 (GPa) 390.0 390.0 401.9
C12 (GPa) 142.7 142 136.4
85
C44 (GPa) 191.0 256 255.7
Young’s Modulus (GPa) 313.6 314.2, 392-694 456.6
Shear Modulus (GPa) 123.7 124, 192 196.5
Poison ratio (𝜈 ) 0.268 0.267, 0.168 0.1616
Table 4.4: Comparison of mechanical properties of 3C-SiC predicted by Vashishta potential with
the experimental and DFT data
Table 4.5 shows the parameters of the Vashishta potential for SiO2 that correctly describes
the structural and mechanical properties of amorphous silica and melting behavior.
Si O
Zi (e) 1.20 -0.60
r1s = 4.43 Å r4s = 2.50 Å rc = 5.5 Å e = 1.60210
-19
C
Two body
Si-Si Si-O O-O
ij 11 9 7
Hij (eV Å

) 0.39246 78.3143 355.5263
Dij (eV Å
4
) 0 3.456 1.728
Wij (eV Å
6
) 0 0 0
Three body
Bijk (eV) o (deg) Cijk  (Å) r0 (Å)
O-Si-O 4.993 109.47 0.0 1.0 2.60
Si-O-Si 19.972 1041.0 0.0 1.0 2.60
Table 4.5: Parameters of the Vashishta interaction potential for SiO2.
86
For the interaction between Al2O3- SiC, SiC- SiO2 and SiO2-Al2O3, we have considered
only the two body terms. These interactions are taken as the geometric mean of the parameters of
the two constituent atoms.
4.2.3 Simulation setup
The initial setup of the RMD simulation is shown schematically in Figure 4.5. It consists
of an alumina matrix with a pre-crack of length 15nm and rows of SiC/SiO2 nanoparticles (NPs)
embedded in the alumina matrix ahead of the pre-crack front. The diameter of these nanoparticles
is between 4nm-7nm with inter-particle distance between 2nm to 5nm. Four sets of simulations
were done corresponding to 1, 2, 3, and 4 rows of NPs. The initial separation of the first, second,
third, and fourth row of NPs from the pre-crack were 4.5 nm, 50 nm, 85 nm, and 115 nm
respectively. Each row had either two or three nanoparticles along the crack front, i.e., the z
direction. The center-to-center distances between nearest neighbor n-SiC was between 2-4 nm
along z and ~ 30 nm along x. The total number of atoms ranged between 6 to 15 million, and the
system size between 40 𝑛𝑚 ×40 𝑛𝑚 ×20 𝑛𝑚 to 180 𝑛𝑚 ×40 𝑛𝑚 ×20 𝑛𝑚 . Equations of
motion in the RMD simulations were integrated with the velocity-Verlet algorithm using a time
step of 1 fs. Systems were first heated to 1,426 C from 25 C in the NPT ensemble using
Berendsen thermostat and barostat in increment of 300 C, where at each temperature it is relaxed
for 20 ps. At 1,426 C it is further relaxed in the NVE ensemble for 100 ps. This process creates
the equilibrated system for the fracture simulations, and periodic boundary conditions were
imposed only during this equilibration process. After equilibration, uniaxial strain was applied in
the system incrementally (0.25%) along the y direction. At each incremental strain, the system was
expanded in the y direction by 0.1 nm over 10 ps, which is then followed by a relaxation step for
60 ps using NVE ensemble. This corresponds to an effective strain rate of 3.5×10
7
sec
−1
.
87
During relaxation phase of the uniaxial compression, Al and O atoms within 2 nm from the
boundaries in the y direction were kept fixed. Under the tensile loading, the pre-notched crack
takes off and propagates along the x axis. In order to avoid spurious effects arising from the
interaction between the crack front and stress waves reflected from the boundaries normal to the x
axis, atoms within 2 nm from the boundaries were held fixed and atoms between 2 – 4 nm from
the boundaries normal to the x axis were damped with Langevin dynamics. A schematic of the
boundary conditions in different regions of the system is shown Figure 4.5b.

Figure 4.5: (a) Schamatic of mode 1 fracture of Al2O3 matrix (light blue) containing a row of n-
SiC (green) with SiO2 shells (red) with a pre-crack (dark blue). (b) Shows the boundary conditions
in different regions of the system.
4.3 Results
4.3.1 Stress-strain curve
We observe that the interparticle distance is very crucial for the crack healing phenomena
as smaller interparticle distance leads to crack healing while large distance between particles is not
effective for crack arrest, and these systems shows brittle fracture like pure Al2O3. In particular,
we observe that the inter-particle distance less than 0.2 nm leads to crack healing. Figure 4.6a
shows the effect of inter-particle distance on stress-strain curve as a function of applied strain for
88
systems consisting of single-row of nanoparticle of size 0.7 nm, and Figure 4.6b shows the effect
of multiple rows on the crack healing and mechanical property of the systems.

Figure 4.6: (a) Effect of inter-particle distance on crack healing Al2O3 composite with single row
of of SiC/SiO2 nanoparticles. (b) Stress-Strain curve for Al2O3 composite with four rows of
SiC/SiO2 nanoparticles and for the pure Al2O3 system during mode 1 fracture. Pure Al2O3 and
systems with larger inter-particle distance fractures at 3.5% however, alumina composite with
nanoparticles withstand up to 5% strain.
From Figure 4.6a-b, we observe that the curves are linear and coincident up to 2% strain
for systems with or without the nanoparticles (NP), however they deviate at higher strain. The MD
results indicate that pure Al2O3 and systems with larger inter-particle distance undergoes brittle
fracture at around 3.5% applied strain. Here, the materials fracture due to the void nucleation and
their coalescence with the advancing crack front in the system. In contrast, alumina nanocomposite
with four rows of SiC/SiO2 core/shell NP with smaller inter-particle distance shows crack healing
and they does not fracture until the strain exceeds 5%.
Our results are consistent with experiments of Ando et al., who introduced indentation
cracks on crack-healed samples and performed three-point bending experiment as a function of
temperature from 25 C to 1525 C. They found fracture outside the crack-healed zones in alumina
nanocomposites. In contrast, fracture initiated just below indentation cracks in samples without
89
crack healing treatment. They also found that bending strengths of crack healed Al2O3 composites
were higher than the bending strength of pure Al2O3 (~1000MPa) at room temperature. Ando et
al. conjectured that the crack-healing mechanism involved the formation of silica shells around n-
SiC by oxidation of Si with atmospheric oxygen followed by diffusion of liquid silica into cracks.
The exothermic heat and volume expansion associated with oxidation of n-SiC are believed to play
important roles in accelerating crack healing by diffusion of silica into the damaged zone. They
have also observed that SiC density inside the Al2O3 matrix is important for crack-healing and has
to be greater than 15% by volume. This is correlated with our MD simulation results that exhibit
interparticle distance and hence nanoparticle density plays an important role in crack healing. Our
simulation results provide a direct evidence for the crack-healing mechanism in Al2O3 containing
silica-coated SiC nanoparticles.
4.3.2 Faceting of alumina cavity
Shapes of cavities in Al2O3 change dramatically under strain at 1,426 C. Figure 4.7(a)
shows a snapshot of an atomic configuration of a cavity containing a n-SiC/SiO2 nanoparticle (NP).
The cavity is no longer spherical. It has facets parallel to prismatic (A) (2
̅
110) prismatic (M)
(1
̅
010) planes of crystalline -Al2O3. These results agree with transmission electron microcopy
(TEM) studies of Hockey et al. on faceting of internal cavities in -Al2O3.
70
Analyzing the
equilibrium Wulff shapes of cavities annealed at 1,600 C, they found well-developed facets in
cavities smaller than 100 nm in diameter. Larger cavities did not exhibit faceting even on
experimental time scales. Cavities are an order-of-magnitude smaller in our computer experiments
and therefore nucleation barriers are small enough for faceting transition to manifest on molecular
dynamics (MD) time scales. This is consistent with theoretical calculations, which indicate that
nucleation barriers for faceting decrease rapidly with a decrease in the cavity size.
71

90
Cavities near the crack facet differently from cavities farther away from the crack. This is due to
the stress field of the crack. When the applied strain increases and the crack advances towards
nearby cavities, the facets of cavities closest to the crack front first shrinks in size and then merge
to change the cavity shape from hexagonal to pentagonal; see Figure 4.7b. These faceting
transitions are caused by surface diffusion of Al and O atoms. Our calculations show that surface
diffusion coefficients of Al and O are 4 ×10
−7
and 2 ×10
−7
cm
2
/s, respectively.
We also observe faceting of n-SiC at 1,426 C, but we do not find any correlation with the
facets of alumina cavities in which nanoparticles are embedded. Kaplan et al. have experimentally
observed faceting of submicron size SiC particles in internal cavities of -Al2O3.
72
Their high
resolution TEM studies reveal that SiC particles develop facets parallel
to {101
̅
1},{1
̅
014} and {1
̅
012} planes and that these facets have no orientational relationship with
facets of alumina cavities which are parallel to {101
̅
2},{101
̅
4} and {0001}

Figure 4.7: Faceting of an Al2O3 cavity around SiC/SiO2 nanoparticle at 3.4% strain. Here, Al2O3
contains two rows of SiC/SiO2 nanoparticle. (a) Al2O3 cavity further away from the pre-crack
changes to hexagonal shape while (b) Al2O3 cavity near the pre-crack changes from spherical to
91
pentagonal shape. These facets form along prismatic (A) (2
̅
110) and prismatic (M)
(1
̅
010) planes.
4.3.3 Crack-healing mechanism

Figure 4.8: (a) Crack healing in Al2O3 by flow of molten silica. Here, the alumina matrix is black
and yellow, and red and green spheres are Si, O and C atoms, respectively. The NP diameters is
7 nm and Al2O3 contains 1,2,3 and 4 rows of NPs. Each row contains 2 NPs of diameter 7 nm
seperated by 2 nm. (b) Self-diffusion coefficients of Si and O as a function of the applied strain.
(c) Strain dependence of Δσxx, which is the difference in the stress at the crack tip and NPs.
Crack healing begins with the diffusion of molten silica through the pinched-off region of
alumina cavities (see Figure 4.7(b)). At the onset of crack healing, silica diffuses through the
nanopores that open up between the crack and the nearest alumina cavities. Unlike bulk, molten
silica which consists of SiO4 tetrahedra, silica in the cavities of alumina at 1,426 C consists of
various SiOx fragments (x = 0, 1, 2 and 3). Figure 4.8(a) presents an atomic view of crack healing
in an Al2O3 composite with a single row of NPs in front of the crack. Here black spheres represent
the Al2O3 matrix and yellow, red and green spheres are Si, O and C atoms, respectively. This
snapshot, taken at a strain of 4%, shows crack blunting and the presence of SiO x fragments on
92
crack surfaces. Figure 4.8(b) shows the effect of the applied strain on self-diffusion coefficients of
SiOx in the direction of the crack propagation. Diffusion coefficients of Si and O are liquid-like,
O diffusion being 25% to 50% faster than that of Si. The diffusion coefficient remains flat and the
crack does not propagate for strains less than 2%. Above 2% strain the crack moves towards
nanoparticles and self-diffusion coefficients of SiOx fragments in the direction of approaching
crack increase linearly with an increase in the applied strain; see Figure 4.8(b). In addition, n-SiC
diffuses into the crack front to stop crack growth. Some of the crack-healing features we observe
are similar to those seen in fluorescent imaging of a crack in a silica film evaporated on PMMA
containing polymer coated CdSe/ZnS NPs.
73
The crack healing mechanism in the experiment
involves diffusion of CdSe/ZnS NPs along the interface between silica and PMMA.
To gain further insight about the crack healing mechanism, we have analyzed the effect of
applied strain on the separation between the crack front and NPs, which indicates how NPs arrest
crack growth. Figure 4.9(a) shows that the separation dC-P between the crack front and two NPs in
the path of the advancing crack decreases with the applied strain. Initially, the size of the two NPs
is 7 nm and the separation between the crack tip and NPs dC-P is 4.5 nm. dC-P decreases to 3.5 nm
as the applied strain yy is increased to 2%. With further increase in strain, dC-P decreases
dramatically and vanishes at yy = 4%. This change in dC-P as a function of applied strain correlates
with the diffusion constant of SiOx fragments from the NPs towards the crack front. It shows that
as the length of the crack front starts increasing, during that time, silica starts flowing towards the
crack and stops its growth.
The applied strain also affects the inter NP distance, dP-P, in much the same way as it does
dC-P, the distance between the crack front and NPs (see Figure 4.9(b)). The initial separation
between the two NPs of diameter 7nm is 1.5 nm. dP-P decreases as the strain is raised to 1.5% and
93
vanishes when the two NPs coalesce at a strain of 2.7%. The two NPs form a neck through which
Si and O atoms begin to diffuse. The width of the neck and Si and O diffusion between the NPs
increase with the applied strain.

Figure 4.9: (a) Distance between crack-front and NPs and (b) distance between NPs dp-p in an
Al2O3 matrix containing a row of two nSiC/SiO2 nanoparticle. The initial separation between the
crack and nanoparticles is 4.5 nm, and the diameter of each nanoparticle is 7 nm.
Changes in stress distribution (Δσ
xx
) between the crack front and NPs provide insight into
crack arrest and healing by Si and O diffusion, and it is the main driving force for silica diffusion
towards the advancing crack. Figure 4.8c shows the stress component Δσ
xx
between the crack tip
and NPs as a function of the applied strain. This stress increases slightly from 4 GPa to 5 GPa
when the strain is increased to 2%. Below this value of strain, Si and O atoms are diffusing in the
shell region of NPs. Further increase in strain creates nanoscale cavities between the crack tip and
NPs. These cavities first coalesce among themselves and then with the crack when the strain
reaches 4%. As a result, Δσ
xx
drops and self-diffusion coefficients of SiOx fragments begin to
increase with the applied strain as these fragments first diffuse into nanocavities and then into the
crack and finally stops its growth.
94
4.3.4 Secondary grain nucleation in strained Al2O3

Figure 4.10: Nucleation and growth of multiple secondary grains at the interface of an Al2O3 –
nSiC/SiO2. Here, the Al2O3 matrix contains two rows of nSiC/SiO2 nanoparticles. Secondary
grains form above 3% strain and facets of the Al2O3 cavity act as nucleation sites for these grains.
Particle 1 in the figure is in row 1, which is closer to the crack then particle 2 in row 2.
At strain above 3%, we also observe formation of secondary grains inside the alumina
matrix which further releases the strain inside matrix and alleviates formation of new cracks inside
the matrix. Figure 4.10(a)-(d) show nucleation of grains (blue) on facets normal to the direction of
crack propagation. These grains are nucleated by the rotation of n-SiC inside the alumina cavities
where facets planes of the alumina cavity as a nucleation site.
74,75
Depending upon whether the
alumina cavity has transformed into a hexagonal or pentagonal shape and also on their location
with respect to the crack front, we observe nucleation of either six or three secondary grains near
the interface of alumina cavity and the nanoparticles. Inside the grains are under coordinated Al
95
(5-fold) and O (3-fold) atoms due to broken Al-O bonds. Grains between neighboring n-SiC grow
with an increase in the applied strains and coalesce at a strain of 4.4%. Grain boundaries are
amorphous, and they contain nanoscale pores at small strains which grow and coalesce to form
secondary cracks at a strain of 4.8%. These cracks are also healed by silica diffusion from cavities.

Figure 4.11: Analysis of stress distribution and local damage inside the Al2O3 matrix containing
2 rows of NPs at 3.4% strain. (a) Stress distribution is plotted on a plane at z=20 nm. (b) Local
damage on a cross-section of Al2O3 at z=20 nm. Here, yellow and red are silicon and oxygen
atoms, respectively. Al2O3 atoms with D
2
min value less than 0.25 are shown in black and the
remaining Al2O3 atoms are colored by their D
2
min value.
We have examined stress distribution in the system to gain further insight into grain
growth. The top panel in Figure 4.11(a) shows the stress component yy in the x-z plane of the
lower half of the system at an applied strain of 3.4%. (yy is time-averaged virial stress in voxels
of size 0.5 nm.
76
) The stress is negligible around n-SiC (see the blue region in the bottom panel)
96
and reaches 15 GPa approximately 10 nm ahead of the first row of n-SiC. The stress decreases
away from the first row of nanoparticles and rises again reaching 12 GPa approximately 10 nm
ahead of the second row of NPs.
The stress yy correlates very well with damage evolution during grain nucleation and grain
growth. Local strain
77
is assessed by calculating local deformation parameter D
2
min for each atom.
D
2
min for each atom is computed by first defining an interaction range (rc) for each atom, which
for alumina is taken as 0.6 nm. After that local strain around each atom is calculated by minimizing
the mean square difference between the actual displacements between the center atom and its
neighbors inside this interaction range, and the relative displacements that they would have if they
were under uniform strain ϵij in a reference frame. Here the equilibrated system at 1,426 C before
the application of tensile force on the system is taken as the reference frame. This minimum value
of ϵij is then used to calculate D
2
min for the central atom i according the equation 1 and the minimum
value of ϵij is computed using equation 2-4.
𝐷 2
(𝑡 ,∆𝑡 )=∑ ∑ {𝑟 𝑛 𝑖 (𝑡 )−𝑟 0
𝑖 (𝑡 )−∑ [(𝛿 𝑖𝑗
+𝜖 𝑖𝑗
)∗(𝑟 𝑛 𝑖 (𝑡 −∆𝑡 )−𝑟 0
𝑖 (𝑡 −∆𝑡 ))]
3
𝑗 =1
}
2
3
𝑖 =1
𝑁 𝑛 =1
Eq 4. 7
where
𝜖 𝑖𝑗
=∑ 𝑋 𝑖𝑘
𝑌 𝑖𝑘
−1
−𝛿 𝑖𝑗 𝑘 Eq 4. 8
and
𝑋 𝑖𝑗
=∑ [𝑟 𝑛 𝑖 (𝑡 )−𝑟 0
𝑖 (𝑡 )]
𝑛 ∗[𝑟 𝑛 𝑗 (𝑡 −∆𝑡 )−𝑟 0
𝑗 (𝑡 −∆𝑡 )] Eq4. 9
𝑌 𝑖𝑗
=∑ [𝑟 𝑛 𝑗 (𝑡 )−𝑟 0
𝑗 (𝑡 )]
𝑛 ∗[𝑟 𝑛 𝑗 (𝑡 −∆𝑡 )−𝑟 0
𝑗 (𝑡 −∆𝑡 )] Eq4. 10
Here r0 and rn are atomic coordinates in the unstrained and strained systems, respectively,
and the summation over n includes all atoms within the range of interaction. D
2
min gives us the
local deviation from affine deformation and can be used as measure of plastic deformation in the
97
system. Its value in unstrained regions is close to zero and regions with large deformation have
value greater than unity.
The bottom panel of Figure 4.11(b) shows the minimized values of Eq. 4.7 along the basal
(0001) and prism (2
̅
110) planes of a section of Al2O3. Here green and blue regions correspond to
D
2
min = 1.5 and D
2
min = 0.25, respectively. These deformation values correlate very well with
changes in yy in the x-z plane. As the increase in size of the regions with high local strain increases
with the applied strain, the tensile stress in those deformed regions drops to zero. We can also
observe from figure 4.11b that the maximum deformation occurs along the slip direction [011
̅
1]
of the prism plane of -Al2O3.
78-80

Structural correlation like the radial distribution functions, g(r), coordination number n(r)
and bond angle distribution shows the difference between the local geometry around each atom in
these secondary grains and pure alumina. Figure 4.12(a) shows two regions of the alumina matrix
in which we have calculated g(r) for Al-O bond, and bond angle distributions for O-Al-O and Al-
O-Al at different values of the applied strain, yy. The formula for the calculation of g(r) and bond
angle distributions is given in equation XX and YY respectively. Figures 4.12(b) and (c) show that
g(r) is nearly the same at yy = 0 and yy = 4.5% in both the secondary grains and pure alumina
however the coordination number of Al (5) and O (3) atoms has decreased slightly in the secondary
grains. In case of bond angle distribution, we do not observe any noticeable difference between
the O-Al-O bond angle distributions at yy = 0 and yy = 4.5% in these two structures. However,
we observe in the Al-O-Al bond angle distribution that the two distinct peaks at 90˚ and 130˚ at
yy = 0 have merged into one broad peak at yy = 4.5% in the secondary grains in comparison of
pure alumina.

98

Figure 4.12: (a) Local deformation analysis of the deformed region in Al2O3 matrix at 4.5% strain.
(b) Comparison of Al-O-Al and O-Al-O bond angle distributions in the deformed region of
alumina. (c) Comparison of radial distribution function g(r) for Al-O and coordination number
n(r) in the deformed region. Here, dotted lines are g(r) and n(r) plots in Al2O3 matrix at 0% strain
and solid lines in the deformed region of Al2O3 matrix at 4.5% strain.
4.3.5 Effect of system size on crack healing
We have also investigated the effect of up to four rows of NPs on damage mitigation in the
alumina matrix. The diameter of each NP (7 nm) and the distance between NPs (2 nm) are kept
the same. Recall, with a single row of NPs, silica heals the pre-crack, but damage appears
approximately 65 nm ahead of the NPs above a strain of 3.5% and alumina matrix fractures at a
strain of 4.5%. Note that from the local deformation analysis we observe that maximum
deformation in the nucleated secondary grains in alumina happens along the [011
̅
1] slip direction
on the prism plane. Because of this, crack develops and propagates on the prism plane along the
slip direction [011
̅
1] due to the formation and coalescence of nanovoids along the [011
̅
1]
direction. In the case of two rows of NPs, cracks appear and propagate directly ahead of the 2nd
99
row of nanoparticles, see Figure 4.13(b). However, the damage zone starts around 100 nm ahead
of the pre-crack. Another row of NPs shifts the damage zone further to 130 nm, see Figure 4.13(c).
Finally, with four rows of nanoparticle the system is healed and there are no cracks in the alumina
matrix at 5% strain, see Figure 4.13(d). The damage is only in the form of nanovoids in between
the rows of n-SiC, but silica diffusion heals them and prevents crack growth. When the strain
exceeds 5%, secondary cracks appear in grain boundaries, but they are also healed by silica
diffusion from nearby cavities.

Figure 4.13: Effect of (a) 2 rows and (b) 4 rows of NPs on crack healing in Al2O3. Here, Al2O3
atoms are shown in black color, and Si and O atoms of NPs are shown in yellow and red,
respectively. Atoms near the crack surface or damaged zone are colored by their D
2
min value.
100
Crack formation happens due to void nucleation and coalescence along the slip direction [01
̅
11]
of the prism plane of -Al2O3, as shown in Figure 4.12(a).
4.3.6 Scalable algorithm to determine fracture surface
During fracture simulation, the determination of fracture surface is important as it gives us
information about the weakest plane in the system and also helps us to make better visualization
of the simulation data. One such algorithm to determine fracture surface is based on connected
component analysis via depth first search (DFS). In this algorithm, we partition the system into
voxels of equal size and based on whether its empty or has atoms, we give it tag. All voxel also
keeps an adjacency list of all its neighbor voxel, and this a graph structure. DFS search on this
graph finds all the voxels which are empty, and this is our fracture surface. But this have large
memory requirement which might become a bottleneck for large simulation data.
By looking at the geometry of the fracture surface, we observe that for each voxel, if we
look at only its neighbors just above and below it, we will be able to determine the fracture. Voxels
near the fracture surface have an empty and a voxel filled with atoms as its neighbors just above
and below it. The only memory requirement this algorithm has the size of the adjacency list for
each node (voxel). Also, identification of each voxels as a part of the fracture is independent from
each other once we compute their neighbor list and hence can computed in parallel using pthreads
and OpenMPI.
4.3.7 Evaluation of mechanical property of Al2O3 after crack-healing by nanoindentation
From the stress-strain curve of alumina and alumina composite with nanoparticles, we
observe that system with nanoparticles shows crack-healing and can withstand strain up to 5%
strain in comparison of pure Al2O3 that fractures at about 3.5% strain. This shows that toughness
of these composite material is higher that of pure Al2O3. However, the strength of these composite
101
systems has slightly decreased. Pure Al2O3 is able to withstand load up to 7.5 GPa whereas crack
healed composites can only withstand load up to 5.0 GPa.
To further gain insight into the mechanical property of Al2O3 after crack healing, we have
performed nano-indentation test of Al2O3 after crack healing at room temperature using SiC
indenter. To perform the indentation test, a region of alumina which is between the two rows of
SiO2 coated SiC nanoparticle in the composite system is taken as initial configuration. The size of
this system was 15𝑛𝑚 ×15𝑛𝑚 ×15𝑛𝑚 , and it is quenched to 25 C under NPT ensemble in 200
ps. After that, it is relaxed at 25 C for another 100 ps. A conical shaped SiC indenter with the tip
angle of 30 is used as an indenter. The indenter is first relaxed using conjugate gradient and after
that heated to 25C under NVT in 50 ps. The indenter is further relaxed for 50 ps at 25 C under
NVE ensemble. After that the indenter is placed at top of the Al2O3 system, and it is moved at a
constant speed of 10 m/sec towards the material. Each indent phase consisted of two stages where
are the displacement of the indenter by 0.1nm in 10 ns which is then followed by the relaxation
phase of entire system for 40 nm where the indenter is held fixed. The total force on the indenter
is computed during the relaxation phase by taking averages over 10 ns. The indenter is moved by
3 nm inside the material, which is then followed by the unloading phase.
To compare the change in mechanical property of the crack healed Al2O3, we have also
performed the indentation on pure Al2O3 crystal. In the two systems, the loading direction was
along the [011
̅
0] direction of Al2O3. Figure 4.14 shows the load-displacement curve for these
systems at room temperature. We c observe that after crack healing the load bearing capacity of
Al2O3 is slightly decreased.
102

Figure 4.14: Load-displacement curve for the nanoindentation of Al2O3 before and after crack
healing at room temperature.
4.4 Conclusion
In summary, RMD simulations of Al2O3 containing SiC/SiO2 nanoparticles show that silica
is an excellent healing agent of cracks in alumina at high temperature. When the applied strain is
increased, bonds break and nanovoids are formed in the alumina matrix between the crack front
and nanoparticles. These nanovoids are the conduits through which silica diffuses into the crack
to stop its growth. The self-diffusion coefficient of silica is liquid-like (> 10
-5
cm
2
/s), and it
increases rapidly with an increase in the applied strain. As a result, the higher the applied strain,
the more rapid the crack healing. RMD simulations also show that alumina cavities containing
SiC/SiO2 nanoparticle are faceted, and the low energy prismatic (A) (2
̅
110) and prismatic (M)
(1
̅
010) planes are the prominent facets. These results are consistent with TEM studies of
nanocavities in -Al2O3. We observe grain formation on facets and grain growth with an increase
in the strain. Damage in the form of nanovoids and nanocracks is observed in grain boundaries,
but rapid diffusion of silica from nanoparticles again heals the damage. The weak link in the
strained Al2O3 matrix is the [011
̅
1] slip direction in the prism plane, along which nanovoids can
103
nucleate and coalesce to form cracks and cause fracture in the [011
̅
1] direction. Nanovoid growth
can be mitigated by a judicious choice of the size and spatial distribution of SiC/SiO2 nanoparticles
along the slip direction.

104
Chapter 5: Strain Induced Phase Transformation in Layered
Materials
5.1 Background
Since the advent of single layer graphene, a wide variety of two-dimensional (2D)
materials, especially semiconducting transition metal dichalcogenides (TMDCs), have been
isolated with a suite of interesting properties and applications.
22-27
These TMDCs have a range of
band gap, for example, the semiconductor with large band gap (MoS2~ 2 eV) and superconducting
PDTe2. Both experiments and simulations have confirmed that the physical properties of these 2D-
materials are sensitive to their layer thickness and composition. This happens due to the spatial
confinement of atoms and charge carriers to atomic thickness in these 2D-materials, and because
of which they exhibit different physical, chemical, electronic and mechanical properties as
compared to their bulk counterpart. Further, we can vertically stack heterogeneous TMDC
monolayers, which opens the door to a whole new class of hybrid materials where one can, in
principle, engineer electronic, transport, optical and other properties in well-defined controllable
material structures.
81-83
An important aspect is the discrete nature of the van der Waals (vdW)
bonded layers that make up these materials combined with the fact that the interlayer interactions
are nonetheless sizable enough to strongly affect electronic behavior. This combination enables
engineering of properties by layer stacking. Hence, they are worthy candidate for next-generation,
functionalized devices due to their tunable band gap, transport and other properties, combined with
mechanical strength and chemical stability. Application of these TMDCs includes as catalyst in
dehydrosulfurization, dry lubricants, photovoltaic power generation and as an ingredient in some
Li-ion batteries.
105
TMDC monolayers have interesting mechanical properties too. For example, Young’s
modulus of MoS2 and WS2 monolayer calculated from the force-displacement curve of
nanoindentation experiment is shown to be ~270 GPa which is higher than the Young’s modulus
of bulk MoS2 (~240 GPa) and WS2 (~170 GPa).
84
Also, several efforts have been made recently
to tune the optical bandgap of the atomically thin sheets of TMDCs by mechanical deformation.
85

Similarly, the alloyed system of the TMCD monolayers also have interesting mechanical
properties. Examples of such materials includes the simultaneous presence of Mo and W
chalcogenides in a monolayer. However, the strain effects on the monolayer TMDC alloys have
not been explored well due to the experimental challenges related to the high sensitivity of TMDCs
under electron beam and limited deformation in atomic force microscopy (AFM).
Experimental Synthesis and mechanical deformation of 2D monolayers
The synthesis of these 2D monolayers are done using mechanical exfoliation, chemical
exfoliation and chemical vapor decomposition (CVD).
86-88
In mechanical exfoliation, thin layer of
TMDCs are peeled from the bulk using scotch tape. After that the cleaved crystal of TMDC on
scotch tape is bought in contact with the substrate (SiO2/Si). Finally, scotch tape is removed, and
we get these thin sheets of TMDCs on the substrate. Ambrosi et at. have chemically exfoliated
MoS2, MoSe2, WS2 and WSe2 using Li intercalation method.
87
The commonly used intercalant for
chemical exfoliation includes n-butyllithium (n-BuLi) and tera-butyllithium (t-BuLi). Whereas in
CVD technique, monolayer samples are grown on an atomically flat single crystal of sapphire
(001) substrate or silica substrate. For example, synthesis of MoS2 using CVD consists of
following steps. Initially powers of MoO3 is placed inside a ceramic boat with the substrate (SiO2,
Si) facing upside down to this ceramic boat. In a separate ceramic boat, sulphur power is placed.
The entire system is placed inside a tube and heated to 650° C in N2 atmosphere. At this
106
temperature, MoO3 is first reduced by S into MoO3-x, which is finally reduced by S into MoS2 as
they diffuse towards the substrate. Figure 5.1a shows the schematics of this MoS2 synthesis
process.
86

Figure 5.1: Synthesis of (a) MoS2 and (b) MoWSe2 monolayer by CVD.
86,89

Similarly, for the synthesis of MoWSe2 monolayer using CVD, a mixture of MoO3 and
WO3 powders is first loaded into a porcelain crucible covered with the sapphire substrates inside
a quartz tube.
88
Another porcelain boat filled with selenium powder is placed upstream next to this
MoO3 and WO3 containing porcelain boat. The tube is filled with a mixture of Ar/H2 gas and
heated to 875 °C for 15 min and cooled down to ambient temperature. The monolayer films are
obtained on the bottom surface of the sapphire. Figure 5.2b shows the schematics of this process.
Mechanical deformation of MoWSe2 monolayer
Mechanical properties (plastic deformation, crack propagation, strain rate effect) of these
TMDC alloys have not been extensively studied due to the experimental challenges associated
with it, like for example high sensitivity of TMDCs under electron beam and their very limited
deformation in AFM. The common approach to study the mechanical properties of TMDC is done
by nanoindentation technique.
89
In nanoindentation, large force and displacement are applied to
the system. The response of the system is then analyzed using the Raman spectroscopy, where we
107
look at the changes in the Raman modes arising from in-plane and out-of-plane oscillatory motion
of atoms due to the applied force. For example, Apte et at. has done the straining experiment of
MoWSe2 monolayer grown on the sapphire substrate.
89
In their experiment, MoWSe2 monolayer
film was first transferred from the sapphire substrate on to a soft PMMA film. After that,
mechanical straining of the MoWSe2 monolayer is done using a nano-indenter. To do that, a three-
point bending method is used to apply the strain and during the experiment the two ends of the
monolayer is firmly clamped. The force on the monolayer is applied using a 2 mm sapphire ball
from underneath using the indentation device. The indenter displacement is converted to strain
values via the following equation, 𝜖 =
𝜏 sin𝜃 2𝑎 . Figure 5.2 shows the schematic of the straining
experiment.

Figure 5.2: (a) Schematic of the mechanical straining of MoWSe2 monolayer by three-point
bending experiment. (b) Shift in Raman mode of MoSe2 and WSe2 as a function of applied strain.
To examine crack propagation in these monolayers, either a sharp or blunt crack is created
in the material first using Berkovich diamond indenter. After that, crack growth is monitored as a
function of applied strain, where strain is again applied using three-point bending test.
Experimental analysis of physical properties using Raman spectroscopy
108
The sensitivity of the physical properties of these TMDC material on layer thickness is
analyzed using Raman spectroscopy (RS). Even the effect mechanical straining and the
composition of the TMDC monolayers is analyzed using Raman spectroscopy. The analysis of
plastic deformation and fracture is simpler with RS as it does not insert defect into the sample that
can make the results noisy. Raman spectroscopy is a spectroscopic technique that is used to analyze
the rotational, vibrational and other low-frequency modes of the materials. It is based on the
principal of inelastic scattering (Raman scattering) of monochromatic light (laser light). Here, the
incident laser light interacts with the molecular vibrations, phonon and other excitation present in
the materials, which shifts the energy of the scattered light. This different in energy between the
incident light and the scattered light is used to identify the vibration modes of the material.

Figure 5.3: Schematic of the atomic displacement of four active Raman-modes in MoS2.
90

The four-active Raman modes observed in bulk TMDCs are E
1
2g, E
1
g, E
1
2g and A
1
g.
90
The
schematics of these modes are shown in figure 5.3. Here E
1
2g, E
1
g, E
1
2g are in-plane mode, and A
1
g.
is out of plane mode of vibration. Out of these four modes, only E
1
g and A
1
g mode is used in the
experimental analysis of TMDCs. Lee et al. have used RS to analyze the thickness dependent
change in frequency shift of Raman modes in MoS2.
90
They have observed E
1
g and A
1
g mode near
384 cm
-1
and 405 cm
-1
respectively. In their experiment, they have observed that in-plane E
1
g
vibration softens (red shift, decreases), and the A
1
g vibration stiffens (blue shifts, increases) with
increasing sample thickness. According to classical model for coupled harmonic oscillator, both
109
E
1
g and A
1
g are expected to increase as a function of layer thickness since the inter-layer Van-der
Waals (VB) attraction increases the net restoring force on the atoms. The unexpected behavior of
E
1
g in their experiment suggests the presence of long-range coulomb interaction or stacking
induced change in interlayer bonding. Not only the frequencies of the E
1
g and A
1
g modes but also
the intensities and the line widths of these two modes show thickness dependence behavior in
MoS2. They have observed that intensities of both of these modes increases linearly with layer
thickness for up to 6 layers, and after that it starts decreasing. In case of line width, it increases
first and then decreases for A
1
g modes; however, E
1
g does not show any thickness dependence.
Selenide TMDCs shows weak in-plane mode (E
1
g) as compared to out of plane mode (A
1
g),
which is different from sulfide where in-plane mode is stronger. In case of MoWSe2 monolayer,
experimentally two peaks are observed at 239 cm
-1
and 249 cm
-1
in RS, which corresponds to the
A
1
g vibration mode of the out-of-plane Se-Mo-Se and Se-W-Se bonds respectively. A weak signal
of E
1
g due to the in-plane vibration mode of Mo-Se-Mo bond is also observed in RS at 285 cm
-1
.
During the experimental synthesis of these alloyed systems, the presence of these peaks in the RS
and X-ray phonon spectroscopy is used to identify the existence of Mo, W, Se in the material.
Raman spectroscopy is also used for the analysis of the mechanical straining of the alloyed system.
As bending strain is applied to the MoWSe2 monolayer, the out-of-plane vibrational modes are
expected to show variations in frequencies, intensities, widths, etc. This happens because the
Raman spectra of the selenide monolayer have stronger out-of-plane vibration mode of Se-Mo-Se
and Se-W-Se bonds as compared to in-plane vibration mode.
The fwhm of a Raman mode indicates the level of plastic deformation and the
instantaneous changes in crystallinity. Analysis of fwhm of the MoSe2 A
1
g and WSe2 A
1
g shows
that Se-Mo-Se bond bears most of the load and hence the fwhm of MoSe2 A
1
g mode is affected
110
more severely than that of WSe2. Due to this mismatched local deformation between MoSe2 and
WSe2 region, a tensile strain acts in the WSe2 region and affects the Se-W-Se bonds. Again, the
local Raman measurements during the crack propagation through the MoWSe2 monolayer shows
that crack propagates through the MoSe2 rich regions, and during crack propagation irreversible
2H to 1T’ phase transformation occurs near the crack tip. This transformation in the experiment is
identified in the monolayer by high-angle annular dark-filed images.
Current Work
Due to the experimental challenges associated with the mechanical deformation of these
TMDC alloys, the mechanism of plastic deformation and crack propagation in these materials have
not be fully explored. Here we have used molecular dynamics simulation to understand the
mechanism of crack propagation and its directional dependence on the loading direction in
MoWSe2 monolayer. To analyze these fracture simulation results, a neural network model is built,
which automatically learns the relevant information of the simulation data like phase
transformation, defect formation.
5.2 Simulation details
5.2.1 Crystal structure of MoSe2 and WSe2
Both MoSe2 and WSe2 are transition metal dichalcogenides (TMDC). Their unit contains
of a single transition element (Mo or W) and two chalcogen atoms (Se). They have hexagonal
crystal structure with space group P63/mmc (194). The unit cell dimensions of MoSe2 and WSe2
are 0.3288 𝑛𝑚 ×0.3288 𝑛𝑚 ×1.2931 𝑛𝑚 and 0.3297 𝑛𝑚 ×0.3297 𝑛𝑚 ×1.2982 𝑛𝑚
respectively. The fractional coordinate of Se and Mo/W in the unit cell of these crystals are
(
1
3
,
2
3
,0.625) and (
1
3
,
2
3
,
1
4
), respectively. These TMDC’s exists in two different polymorphs, which
are H and T’ crystal structure. In H crystal structure, each transition metal has trigonal prismatic
111
arrangement of chalcogen atoms around it and lacks inversion symmetry. These H polytypes are
direct band-gap semiconductor. Whereas, in T’ crystal structure, each transition metal has either
distorted or twisted octahedral symmetry of chalcogen atoms around it and have inversion
symmetry. The T’ polytypes have semi-metallic electronic structure. However, in both of these
crystal structure, each transition element is surrounded by 6 chalcogen atoms, and each chalcogen
atoms is surrounded by 3 transition elements in a monolayer.
Controllable transformation from 2H to 1T’ phase is of central importance to electronic
industries as it will help in the formation of defect-free and mechanically robust semiconductor-
metal hetero-phase junctions. Currently the common techniques for this 2H to 1T’ transformation
are based on thermal annealing and mechanical straining.
91

5.2.2 Stillinger-Weber interaction potential for MoWSe2
Stillinger-Weber (SW) is a semi-empirical force field originally proposed by Stillinger and
Weber to describe the interaction in solid and liquid forms of silicon.
92
After that it has be
developed for many other covalent materials like MoS2, WSe2 and single layer black phosphorous.
It has been extensively used in MD simulations for thermal and elastic phenomena of transition
metal dichalcogenide (TMDC) materials. Its functional form consists of a 2-body and a 3-body
term and is written as
𝐸 =∑ 𝑉 2
(𝑟 𝑖𝑗
)
𝑁 𝑖 +∑ 𝑉 3
(𝑟 𝑖𝑗
,𝑟 𝑗𝑘
,𝜃 𝑖𝑗𝑘 )
𝑁 𝑖 Eq5.1
𝑉 2
(𝑟 𝑖𝑗
)=𝐴 ⋅exp(
𝜌 𝑟 𝑖𝑗 −𝑟 𝑐𝑢𝑡 )(
𝐵 𝑟 𝑖𝑗 4
−1) Eq5.2
𝐸 𝑎𝑛𝑔𝑙𝑒 (𝑟 𝑖𝑗
,𝑟 𝑗𝑘
,𝜃 𝑖𝑗𝑘 )=𝜆 ⋅exp(
𝛾 1
𝑟 𝑖𝑗 −𝑟 𝑐𝑢𝑡 +
𝛾 2
𝑟 𝑗𝑘
−𝑟 𝑐𝑢𝑡 )(cos𝜃 𝑖𝑗𝑘 −cos𝜃 0
)
2
Eq5.3

112
Here, E is the total energy of the system, and V2 and V3 are contribution due to the bonding
and the angular terms respectively to the total energy. rij and rjk denotes the bond length between
atom i and j, and j and k, respectively. 𝜃 𝑖𝑗𝑘
is the angle between the atoms i, j and k, where j is the
central atom, and 𝜃 0
is the equilibrium angle. The exponential factors in the two-body and three-
body terms results in smooth and rapid decay of the interaction potential, and hence force and
energy computations are very fast in SW potential. A, B, 𝜌 , 𝛾 1
, 𝛾 2
and 𝑟 𝑐𝑢𝑡 are tunable parameters
that are fitted to get the correct lattice constant, cohesive energy, phonon dispersion curves and
elastic constants. Here 𝑟 𝑐𝑢𝑡 is the cutoff distance of the potential.
For the MoSe2 system, we use the Stillinger Weber force field generated by Kandemir et.
al., which was parametrized against elastic and structural properties of the MoSe2 crystal.
93
The
force field parameters for the WSe2 crystal is optimized using density functional theory
calculations of structural properties like lattice constants and mechanical properties like Young’s
modulus and shear modulus. Force field parametrization was performed with the GULP code.

Table 5.1: 2 body terms Stillinger Weber parameters for the MoWSe2 crystal.

Bond 𝑨 𝝆 𝑩 𝒓 𝒄𝒖𝒕
Se1-Se1 1.841 0.100 20.000 4.05
Se2-Se2 1.841 0.100 20.000 4.05
Se1-Se2 1.841 0.100 20.000 4.05
Mo-Se1 30.000 2.999 30.000 3.36
Mo-Se2 30.000 2.999 30.000 3.36
Mo-Mo 3.374 1.000 38.267 4.55
W-Se1 47.118 2.999 30.000 3.36
W-Se2 47.118 2.999 30.000 3.36
W-W 3.706 1.000 38.267 4.55
Mo-W 3.536 1.000 38.267 4.55
113
The Stillinger Weber parameters for MoWSe2 crystals used in this study are listed tables
5.1-5.2. For the MoSe2 and WSe2 interface, we have considered the following terms: Mo-W
bonding and thr angular contribution from Mo-Se1-W, Mo-Se2-W. Here, Se1 and Se2 are Se atoms
in top and botton layers respectively.
Angle 𝝀 𝟏 ,𝝀 𝟐 𝜽 𝟎 𝜸 𝒓 𝒄𝒖𝒕
Se1-Mo-Se1 20.000 81.530 2.000 3.36
Se2-Mo-Se2 20.000 81.530 2.000 3.36
Mo-Se1-Mo 20.000 81.530 2.000 3.36
Mo-Se2-Mo 20.000 81.530 2.000 3.36
Se1-Mo-Se2 13.128 82.128 2.000 3.36
Se1-W-Se1 33.838 81.530 2.000 3.36
Se2-W-Se2 33.838 81.530 2.000 3.36
W-Se1-W 17.512 81.530 2.000 3.36
W-Se2-W 17.512 81.530 2.000 3.36
Se1-W-Se2 22.212 82.128 2.000 3.36
Mo-Se1-W 18.715 81.530 2.000 3.36
Mo-Se2-W 18.715 81.530 2.000 3.36
Table 5.2: 3 body terms Stillinger Weber parameters for the MoWSe2 crystal.
Material Quantity SW forcefield Experimental/DFT
MoSe2 Young’s modulus (GPa) 152.5 162.1
Lattice constant (a, Å) 3.31 3.3
Included unit cell angle (γ, °) 120 120
WSe2 Young’s modulus (GPa) 160.9 165

Shear modulus (GPa) 32.0
26.8
39.2
Lattice constant (a, Å) 3.31 3.28
Included unit cell angle (γ, °) 120 120
114
Table 5.3: Comparison of Stillinger Weber properties with DFT and experimental values for
MoSe2
94,95
.
Table 5.3 gives the comparison the structural and elastic properties of MoSe2 and WSe2
calculated by the optimized SW force field to values obtained from experiments and DFT
simulations.
5.2.3 Simulation setup
Molecular dynamics (MD) simulations are performed to examine vibrational modes and
atomistic mechanisms of defect formation and crack propagation in the strained MoWSe2 alloy.
The system consists of an MoSe2 monolayer containing regions of either a single patch of WSe2
or multiple patches of WSe2 distributed randomly in the MoSe2 matrix. For systems with multiple
patches of WSe2, the fraction of WSe2 in MoSe2 is kept same as in the experiment, which is 30%.
Simulations are performed for two different system sizes which are .5 𝜇 m x 0.5 𝜇 m and 1 𝜇 m x 1
𝜇 m. In all systems, the MoSe2 matrix has a pre-crack of length 0.166 𝜇 m.
After creating the systems, they are relaxed using conjugate gradient method and then
heated to 100 K under NVT ensemble in 100 ps. At 100 K, it is further relaxed for 100 ps in NVE
ensemble. After relaxation, it is subjected to homogeneous tensile force. Figure 5.4 shows the
schematic of mode 1 fracture simulation of MoWSe2 monolayer. The loading direction of the
tensile force is either along the zigzag (10) or the armchair (11) direction of the crystal structure,
such that the pre-crack is perpendicular to the loading direction. In (10) direction, applied tensile
load is along the Mo-Se bond whereas in (11) direction it is perpendicular to the Mo-Se bond as
shown in figure 1b-c.
Each deformation phase consists of two stages, which are homogeneous tensile
deformation followed a relaxation phase. During the homogeneous tensile deformation phase, the
system size is expanded homogeneously along the loading direction by 0.225% in 100 ps which is
then followed by the relaxation phase where the entire system is relaxed under NVE ensemble.
115
The whole process is repeated until the system fractures, and the effective strain rate is 2.25  10
7

s
-1
. Also, a 0.2 nm region near the boundaries are clamped for the entire time, and force along the
z-direction is made zero both during the relaxation of the system and during the application of
tensile forces. No periodic boundary condition was applied in any direction.

Figure 5.4: (a) Schematic of mode 1 fracture in MoWSe2 alloy containing a single patch of WSe2.
(b) and (c) shows atomic view when loading direction is perpendicular to (1,1) and (1,0),
respectively.
5.3 Mode-1 fracture of MoWSe2
5.3.1 Fracture toughness of MoSe2 and WSe2
First MD simulations are performed to calculate the fracture toughness K1c of MoSe2 and
WSe2. To do that stresses were calculated around the pre-crack tip and fitted to the results of linear
elastic fracture mechanics outside the plastic deformation zone. According to linear elastic fracture
mechanics, the stress distribution near the crack tip in polar coordinates (𝑟 ,𝜃 ) with origin at the
crack tip follows the following relation
116
𝜎 =
𝐾 √2𝜋𝑟
[cos
𝜃 2
(1+sin
𝜃 2
sin
3𝜃 2
)] Eq5. 4
Which for 𝜃 =0 becomes
𝜎 =
𝐾 √2𝜋𝑟
Eq5. 5

Figure 5.5: KIC value along different loading direction in MoSe2 and WSe2
Here, 𝜎 is normal stress along the loading direction in the material in-front of the crack tip.
r is distance from the crack tip and K is the proportionality of constant known as stress-intensity
factor. The value of K at the critical stress when the system fracture is the fracture toughness (KIC)
of the material. KIC is the materials property and is independent of the crack length, geometry and
loading direction. It describes the ability of the material to resist fracture and is higher for ductile
material as comparted to brittle materials. Figure 5.5 shows the KIC value for MoSe2 and WSe2 for
(1,0) and (1,1) direction. It is computed at the critical strain just before the fracture by plotting the
117
stress in front of the crack tip as a function of distanced for 𝜃 =0 and then fitting equation 5.5
whose slope is the KIC. The computed value of KIC matches with the experimental value of K1c that
were obtained by using Griffith’s criterion in conjunction with measurements of surface energies.
Further, we can observe that WSe2 has higher fracture toughness than MoSe2 and in MoSe2, KIC
is lower along the (1,0) direction. Also, lower KIC value indicates that both of these materials are
brittle in nature as suggested in the literature.
5.3.2 Mode-1 fracture of MoWSe2 along (10) direction

Figure 5.6: Distribution of the stress component yy during crack propagation in MoWSe2
monolayer for (a) (1,0) and (b) (1,1) loading direction. Critical stresses near the crack tips reach
20 GPa at the onset of fracture.
The system fracture at strain above 2% which matches with the experimental observation
(~2.16 % strain). Figure 5.6 shows the stress distribution in the system during crack propagation.
The maximum stress at the crack tip just before the crack propagation reach up to 20GPa. The
regions near the crack tip have damaged areas due to broken Mo-Se bonds. The bonds length of
118
Mo-Se and W-Se is increased to 0.27 nm and 0.268 nm respectively in the strained system, which
is higher than the equilibrium bond length of Mo-Se (0.253 nm) and W-Se (0.254). This means
that Mo-Se bonds are more strained as compared to W-Se regions, which is consistent with the
mechanical straining experiments of MoWSe2 monolayer that shows that Se-Mo-Se bond bears
most of the load. Here the material shows cleavage fracture, and crack propagates directly through
the material. The observed crack speed during crack propagation is 2.00 km s
-1
.
5.3.3 Mode 1 fracture of MoWSe2 along (11) direction
MoWSe2 monolayer fractures at a strain above 2.15 % when the loading direction is along
(1,1). Here the maximum stress at the crack tip just before the crack propagation reach up to 20
GPa. Figure 5.6b shows snapshots of crack propagation and process zone just after the onset of
crack propagation in the system. Unlike (10) loading direction, here the crack meanders slightly
as it propagates through the MoSe2 matrix by breaking Mo-Se bonds. The pre-crack branches into
daughter cracks as it crosses the MoSe2/WSe2 interface and enters the WSe2 patches (Figure 5.6b,
5.8). The most dramatic change occurs in the process zone around the crack tip, where large stress
concentration and the resulting biaxial tensile strain trigger an irreversible local structural phase
transformation (Figure 5.7a) from the ground state 2H crystal structure to the 1T crystal structure.
This is consistent with previous theoretical predictions and experimental observations in other
elastically strained TMDCs
17,18
. This irreversible 2H to 1T phase transformation during crack
propagation not only for the (1,1) loading direction but also for the (1,0) loading direction.
119

Figure 5.7: Stress induced 2H-1T phase transformation near the tip of the propagating crack. (a)
Snapshot of the 2H-1T phase transformation as the crack propagates through the MoSe2/WSe2
alloy. The highlighted region indicates the two different crystal structures in the process zone.
Green atoms are at the interface of the 2H phase and the disordered 1T phase. (b) Simulated
Raman spectra of the MoSe2/WSe2 alloy during crack propagation. Spectra obtained during the
initial stages of crack propagation (0-1 ps) are characteristic of the 2H crystal structure of the
MoSe2/WSe2 alloy, while spectra from later times are obtained from simulation cells with greater
1T content. Figures (c)-(e) are snapshots from the MD simulation of atomic configurations in the
process zone. These snapshots show atomic displacements corresponding to the motion of the 2H-
1T phase transformation front (dashed brown line). Here Mo, Se, and W are shown in pink, yellow
and blue colors, respectively. The transformation of Mo atoms (1, 4, 7) to the 1T structure is due
to the displacement of Se atoms (2, 3, 5, 6, 8 and 9) highlighted in figures (c)-(e). Dotted blue
curve indicates bond breaking and dotted blue lines shows bond formation. For example, in figure
(c) the formation of 1(Mo)-5(Se) bond and breaking of 1(Mo)-2(Se) and 7(Mo)-5(Se) bonds
transforms 1(Mo, cyan) from 2H to 1T structure. The resultant 1T structure is shown in figure (d).
During this transformation, the Se pairs 2-3 and 5-6 move in opposite direction as shown in figure
(d).
120
Figure 5.7b shows the simulated Raman spectra from the 2H and 1T regions of the
MoWSe2 monolayer. Raman spectra are obtained from a calculation of the Fourier transform of
the group breathing autocorrelation functions. These autocorrelation functions are calculated from
the atomic trajectories in the MD simulation of the 0.5 µm ⨉ 0.5 µm monolayer by projecting the
relative velocity between each Mo/W cation and its neighboring Se anion along the direction of
the cation-anion bond. This dynamical property is associated with bond stretching and can be used
to accurately determine the frequencies of the Raman absorption peaks, however spectral
intensities obtained from this method are unreliable, since information about Raman scattering
cross sections are not considered. Time-resolved Raman spectra obtained during crack propagation
through the MoSe2/WSe2 alloy monolayer show that spectra from the initial stages of crack
propagation are characterized by peaks at 220 cm
-1
and 270 cm
-1
. Spectra obtained from after crack
propagation through the sample reflect the stress redistribution and phase transformation through
the emergence of a new intermediate peak at 250 cm
-1
. This trend is also observed in the
experiment where the Raman spectra at the crack tip and edges of the strained alloy shows
convoluted peaks between 250 and 260 cm
-1
.
5.3.4 Mechanism of strain induced 2H to 1T’ phase transformation
Crack propagation in the alloy is accompanied by a structural phase transformation from
2H and 1T phase. These common polymorphs in the TMDC family of materials differ in the
arrangement of chalcogen atoms around the central transition metal atom. The 2H crystal structure
results from a A-B-A-type stacking in the chalcogen-transition metal-chalcogen trilayer leading to
a trigonal prismatic coordination around the central Mo ion (Figure 5.7a, bottom inset). The 1T
crystal structure arises from an A-B-C-type stacking of chalcogen-transition metal-chalcogen
trilayer resulting in a distorted octahedral coordination around the central Mo ion (Figure 5.7a, top
121
inset). The 1T crystal structure, with a larger lattice constant, is metastable, but becomes
energetically favorable under high biaxial tensile strains. This stress relieving 2H-1T phase
transformation proceeds primarily via the in-plane gliding of one of the chalcogen layers in the
armchair direction (Figure 5.7c, e), identical to atomic mechanisms observed in 2H-1T phase
transformations in other TMDC materials
19
. Experimental evidence of stress induced phase
transitions in other ternary TMDC alloys has been observed. For example, effect of applied strain
near holes in MoWS2 alloys showed that regions with structural phase transformation near the
holes becomes larger as holes grow under applied strain.

Figure 5.8: Local deformation around the crack tip in a strained in an MoSe2/WSe2 monolayer in
[1,1] configuration. (a) Local strains in the process zone around the pre-crack and crack
branches. In an unstrained MoWSe2 monolayer in the 2H or 1T phase, Mo and W atoms belong to
six and Se to three rings, each with four members. Here, green dots represent defects, i.e., Mo and
Se atoms in which one of the 4-membered ring is replaced by a larger ring.
To further characterize local strain in the system during crack propagation, we also
calculate strains from MD trajectories using a local deformation analysis technique as shown in
figure 5.8. We can observe that the local strain is higher in nearby region (process zone) of the
122
crack where structural transformation has occurred. This process zone also has atomic and
extended defects, which we have analyzed using the the ring structures of Mo, W and Se atoms.
In an unstrained MoWSe2 monolayer, every Mo, W and Se atom belongs to a fixed number of
four-member ring: 6 rings in the case of Mo or W and 4 in the case of Se. In the process zone
shown by green dots in figure 5.8, many W atoms belong to five and Se to three rings with four
members and a larger ring due to broken W-Se bonds. These defects are distributed throughout the
process zone, more so at the interface between the 2H and 1T phases (see figure 5.8). In addition
to these atomic defects, secondary cracks (red lines) nucleate from the sides of the pre-crack and
they propagate through the 1T phase of the process zone. As the following discussion shows, these
point defects play an important role in controlling the crack branching mechanism and therefore
significantly affect the toughness of the monolayer material.
5.3.5 Mechanism of crack branching and crack closure
As crack propagates through the MoWSe2 monolayer, it branches into multiple daughter
crack for (1,1) loading direction. Some of these branches grow while others crack healing by crack
closure. Hence, the crack-branching mechanism offers insight into the increased toughness of the
MoWSe2 alloy. Crack propagation occurs through the cleavage of Mo-Se and W-Se bonds at the
crack tip due to the intense stress concentration (with σyy approaching 16 GPa at the MoSe
2
/WSe
2

interface, as seen in Figure 5.6b). Figure 5.9a-b shows successive snapshots of atomic
configurations near the crack tip in the WSe2 patch from the MD simulations, which demonstrate
that crack growth leads to the nucleation of point defects in the process zone on either side of the
propagating crack tip. These defects subsequently coalesce to create a preferred path for crack
propagation, leading to the formation of multiple crack branches in the 2D monolayer.
123
The velocity of this process is not uniform in all direction around the crack tip. The average
velocity of the process zone front (defined as the boundary between the 2H and 1T phases) is
between .2 ± 0.6 km s
-1
for these alloyed systems. In case (1,1) loading, the velocity of process
zone along the crack growth direction (perpendicular to the loading direction) is faster than the
crack speed whereas for (1,0) loading reverse is true. Since the nucleation of crack branch happens
inside the process zone and due to this reason crack branching not observed for (1,0) direction.

Figure 5.9: (a-d) shows crack propagation and crack branching mechanism in MoSe2/WSe2
monolayer. (a) At the crack tip, multiple W (orange)-Se (cyan), W (orange)-Se (light blue) and W
(pink)-Se (light blue) bonds breaking leads to the formation of defects at the crack tip along
different directions as shown in (b) which eventually leads into (c-d) crack branching in the
material.
124
5.3.6 Toughing mechanism in MoWSe2 hetero-structure

Figure 5.10: Area of the process zone in the MD simulation as a function of time in a pure MoSe2
monolayer (blue) and in a MoWSe2 alloy (red).
Mode 1 fracture simulation is also done for pure MoSe2 so as to compare the increased
toughness of these alloyed systems. Firstly, the rate of crack growth, as measured from MD
simulations on systems of size 0.15 µm ⨉ 0.15 µm, is smaller in the case of the MoWSe2 alloy
(crack velocity = 1.95 km s
-1
) than in the case of the pure, unalloyed MoSe2 monolayer (crack
velocity = 2.2 km s
-1
). Second, the area of the process zone accompanying the growing crack is
smaller in the MoWSe2 alloyed system as shown in figure 5.10. The average velocity of the process
zone front (defined as the boundary between the 2H and 1T phases) was calculated to be 2.7 ± 0.5
km s
-1
in the pure MoSe2 monolayer, and 2.2 ± 0.6 km s
-1
in the MoWSe2 alloy for (1,1) loading.
Consistent with the anisotropic distribution of local stresses in the process zone, this process zone
growth is much faster along the direction of crack growth. The measured process zone velocity in
front of the crack is 3.5 ± 0.3 km s
-1
in the pure MoSe2 monolayer and 2.6 ± 0.3 km s
-1
in the alloy,
which is consistent with the increased toughness of the MoWSe2 alloy. In addition to this, crack
125
closure, nucleation of secondary grains in the 1T structure of the strained MoSe2 matrix are
additional toughening mechanism in MoSe2 and MoWSe2 monolayers.
5.3.7 Summary
In conclusion, MD simulations of mode 1 fracture in MoWSe2 monolayer reveal that crack
growth in the alloy is accompanied by a structural phase transformation from the 2H to the 1T
phase in the process zone around the crack. Crack propagation mechanism in these systems is
highly directional dependent, where (1,1) loading shows crack branching and crack closure,
whereas in (1,0) loading it shows cleavage fracture. Toughing mechanism during crack
propagation includes structural phase transformation from the 2H to the 1T phase, crack branching,
crack closure, secondary grain nucleation. Also, these alloyed system shows higher toughness as
comparted to pure systems due to lower crack speed and process zone volume
5.4 Neural network analysis of dynamic fracture in a layered material
5.4.1 Introduction
Molecular dynamics (MD) simulation of various physical and chemical phenomena of
materials requires complex data analysis of the MD simulation results so as to identify different
phases, chemical reaction and defects generated during the simulation. For example, during
nanoindentation simulation of metallic systems, plastic deformation happens inside the system due
to dislocation nucleation on (111) plane. Similarly, under stress, materials like SiC, AlN show
phase transformation. Identification of different phases and defects generated during MD
simulation requires complex structural analysis of each atom. These analyses can be as simple as
the calculation of nearest neighbors or the identification of all shortest circuit for each atom. In
general, there are no single evaluation technique that can be applied to all simulation data since
these analyses are very data dependent. However, we can observe that the different structures
126
generated in the MD simulation data are nothing but some complex functional form of the local
environment of the atoms. Hence, we can create a machine learning (ML) model that can
automatically learn this functional directly from the simulation data.
ML methods has been used recently in the material science domain and has shown
phenomenal success for high-throughput screening and property prediction.
96-98
For example,
statistical models built using a smaller fraction of the crystal structures (called training data) can
accurately predict wide range of materials property like band gap,
99,100
dielectric breakdown
strength
101,102
and melting point
100
for remaining crystal structure without doing any quantum
mechanical (QMD) calculation for them. ML methods like support vector regression
99,100
, neural
network
103
and kernel ridge regression
104
are shown to accurately predict the bandgap of double
perovskites, polymers
105,106
and defects in FCC and disordered solid. Neural network and gaussian
processes model (GP) are also used to create data driven force field for materials, where parameters
of the model are learned using the QMD simulation trajectory. These models are highly robust and
can automatically tune their parameters on the fly whenever there are uncertainties.
To identify the different phases generated during the dynamic fracture in MoWSe2
monolayer, we have trained a three-layer, feed-forward neural network (NN) model. The training
data set consists of 36,000 atoms, where each atom is represented by a 50-dimension feature vector
consisting of radial and angular symmetry functions. Hyper parameters of the symmetry functions
and network architecture are tuned to minimize the model complexity with high predictive power
using feature learning and Bayesian optimization. The mode classifies each atom in one of the six
phases which are either as transition metal or chalcogen atoms in 2H phase, 1T phase and defects.
Comparted to tradition approach of structure identification in layered materials, which requires
more than one structural analysis like shortest circuit analysis and bond angle distribution for each
127
atom, NN model is highly scalable and can used to analyze large data set quickly. Further, the
class probabilities can be used to make richer visualization of data with variant colors compared
to traditional method, which assigns single color to each phase.
5.4.2 Model building
To build the NN model, we need to first define a mathematical representation for each
atom. This is usually a high-dimension vector called feature vector that is representative of its local
environment. There are many choices for the feature vector of atoms, which includes
electronegativity (EN), ionization potential (IP), atomic displacement, atomic mass, local stress,
radial and angular symmetry functions. In general, the choice of right feature vector is highly
system specific. For example, EN, IP and atomic mass is suitable choice to build a regression
model for band gap for different class of materials whereas, radial and annular symmetry functions
are suitable to create NN and GP force fields. Here, we have represented each atom by 436-
dimension feature vector, which consists of combinations of radial and angular symmetry
functions. The equations for radial and angular feature are shown in equation 5.6 and 5.8
respectively.
𝑓 𝑐 𝑅 (𝑅 𝑖𝑗
)={
∑ 𝑒 𝜂 1
(𝑅 𝑖𝑗
−𝑅 𝑐𝑢𝑡 )
2
.𝑓 𝑐 (𝑅 𝑖𝑗
) 𝑅 𝑖𝑗
<𝑅 𝑐𝑢𝑡
𝑁 𝑗 =1
0 𝑅 𝑖𝑗
≥𝑅 𝑐𝑢𝑡
Eq5. 6
𝑓 𝑐 𝐴 (𝑅 𝑖𝑗
)={
2
1−𝜁 ∑ (1+𝜆𝑐𝑜𝑠 𝜃 𝑖𝑗𝑘
)
𝜁 .𝑒 𝜂 1
(𝑅 𝑖𝑗
2
+𝑅 𝑖𝑘
2
)
.𝑓 𝑐 (𝑅 𝑖𝑗
).𝑓 𝑐 (𝑅 𝑖𝑘
)
𝑁 𝑗 =1
𝑅 𝑖𝑗
<𝑅 𝑐 𝑢𝑡
0 𝑅 𝑖𝑗
≥𝑅 𝑐𝑢𝑡
Eq5. 7
𝑓 𝑐 (𝑅 𝑖𝑗
)={
0.5 ×[cos(
𝜋 𝑅 𝑖𝑗
𝑅 𝑐 )+1] 𝑅 𝑖𝑗
<𝑅 𝑐𝑢𝑡
0 𝑅 𝑖𝑗
≥ 𝑅 𝑐𝑢𝑡
Eq5. 8
Here in equation 5.6-5.8, 𝜂 1
,𝜂 2
𝜁 ,𝜆 and Rcut are tunable parameters, which is also known
as hyper-parameter of the model. These parameters determine the width and the mean value of the
128
symmetry functions. For each atom i, contribution due to all its N neighbors that are within the
cutoff distance Rcut (8Å) is computed for a specific value of radial (𝜂 1
) and angular (𝜂 2
, 𝜁 , 𝜆 )
features, respectively. Effect of neighbor atoms that are further away from the central atom is
damped using equation 5.7. The computed value of 𝑓 𝑐 𝐴 (𝑅 𝑖𝑗
) and 𝑓 𝑐 𝐴 (𝑅 𝑖𝑗
) for each hyperparameter
set is stored along one of the dimensions of the feature vector. The values of the hyper-parameter
are chosen in each a way that they completely capture the local neighbored information. Here,
radial feature vector for each atom is created due to the contribution from Mo-Mo, Mo-Se and Se-
Se bond where 𝜂 1
is varied between 0.0001 to 1000 on a log scale. For the angular features,
contribution from Mo-Se-Mo, Se-Mo-Se, Mo-Mo-Se, Se-Se-Mo bonds are used. For the angular
features, 𝜆 is taken as -1 or +1, 𝜁 and 𝜂 2
between 0.1 to 10. In total, we have 113 radial feature and
323 angular features for each atom, and each feature has its own tunable hyper-parameter.
However, all features are not equally important. For example, in layered materials, there are no
bond angle of value 50º or distance between two atoms less than 1Å. Instead of removing these
features manually, which is a cumbersome process as all the irrelevant features are not known
beforehand. Thus, features corresponding to the irrelevant attributes are removed using feature
selection, and in the final model each atom is represented by a 50-dimension vector. Further the
hyper-parameters like Rcut for radial and angular symmetry function is optimized using Bayesian
optimization, which increases the model accuracy up to 95%.
After creating the feature vector for each atom, feed-forward neural network with three
hidden layers is trained to identify the different phases of MoWSe2 monolayer. Output layer of the
NN model consists of 6 classes, which are Mo in 2H phase, Mo in 1T phase, Mo as defect, Se in
1H phase, Se in 1T phase and Se as defect. Here, input to NN model is the feature vector of the
atoms and the output is their probabilities in these six classes, and we choses the class with
129
maximum probability as the phase label for the atoms. The first, second and third hidden layer has
350, 100 and 50 hidden units. In the first and second hidden layer, we have used Relu
(y=max (0,𝑥 ) ) activation function while in the third hidden layer Sigmoid (𝜎 (𝑥 )=
(1+𝑒 −𝑥 )
−1
) is used. Figure 5.11 shows a schematic of the NN model, which shows the
conversion of atoms to feature vector and all the six output classes of the NN model. The NN
model is trained using 36,000 examples, which consists of equal distribution of atomic data from
all six classes, and a fracture simulation data of 100,000 atoms is used as a validation set. The
model is trained by minimizing the SoftMax (equation 5.9) function using Adam-optimizer.
Training of NN model is done for a total of 200 epoch, where in each epoch a batch size of 512 is
used with a dropout probability of 0.25 for each hidden layer.
𝜎 (𝑍 )
𝑗 =
𝑒 𝑗 𝑧 ∑ 𝑒 𝑗 𝑧 𝑘
𝐾 𝑘 =1
𝑓𝑜𝑟 𝑗 =1,2,…6. 𝐻𝑒𝑟𝑒 𝑗 𝑎𝑟𝑒 𝑡 ℎ𝑒 𝑜𝑢𝑡𝑝𝑢𝑡 𝑐𝑙𝑎𝑠𝑠 Eq5. 9

Figure 5.11: Schematic of the neural network model showing the mapping of atomic coordinates
into sets of radial and angular feature vectors, NN layer architecture and the six out phases
predicted by the model.
130
5.4.3 Visualization of the learned data representation of NN model

Figure 5.12: (a-c) shows PCA and (d-e) shows t-SNE analysis of the input feature vector and the
learned representation of data in hidden layer 1 and 3 of the neural network.
Figure 5.12 shows the learned representation of the training data set by the first and third
layer of the neural network and the original feature space by principal component analysis (PCA)
and t-distributed stochastic neighbor embedding (t-SNE). Both PCA and t-SNE are dimensionality
reduction technique and they both represent the high-dimensional data in lower dimension.
However, PCA preserves the global information like the distance between different clusters,
nearest neighbor of each clusters in the original high-dimension but in doing so its loses the local
information of the data. For example, here in figure 5.12a-b, we can observe that PCA is able to
distinguish the two major phases which are transition metal (Mo, W) and chalcogen atoms (Se)
but it fails to the captures the local variation in their structural representation, which are 2H crystal,
1T’ crystal and defects. The separation between 2H and 1T’ crystal structure become clearer in the
third hidden layer as it contains four clusters compared to hidden layer 1,2 and the original feature
space who contains only two clusters. Here, the axis of figure 5.12a-c corresponds to the largest
two eigenvalues of the feature vector of the training data set.
131
In comparison of PCA, t-SNE is a non-linear dimensionality reduction technique to
visualize high dimensional data in low dimension. It does not preserve global information of
different clusters, but it is able capture the local variation of the high-dimensional data by
preserving the local neighborhood information of the high-dimensional data in the lower
dimensional representation. It does this by minimizing the KL-divergence between the two
distributions. The expression for this objective function (KL-divergence minimization) is given as
𝐾𝐿 (𝑄 ||𝑃 )=∑ 𝑄 𝑖𝑗
(
𝑄 𝑖𝑗
𝑃 𝑖𝑗
)
𝑖𝑗
Eq5. 10
Here, Pij and Qij are the distribution of the N data points in the higher and the lower
dimensions, respectively. Since KL is an asymmetric function, it gives large cost for representing
nearby data points in higher dimension by widely separated those points in the lower dimension,
and lower cost for representing widely separated points in higher dimension by nearby points in
the lower dimension. Hence, in doing so it is able to preserve the local structure. t-SNE also have
a hyper-parameter called perplexity, which determines the number of nearest neighbors for each
data point, and it needs to be tuned. A low value of perplexity considers every data point as separate
cluster whereas a very high value perplexity considers entire data set as one big cluster. In general,
the perplexity is taken a number between 5 to 50.
Figure 5.12d-f shows the t-SNE representation of training data set in 2-dimension. We can
observe that separation between all the six phases become clearer in the third hidden layer as
compared to in hidden layer 1 and the original feature space. In the original feature vector
representation of atoms, each atom is represented as some combination of radial and angular
symmetry functions, which does not the non-linear interaction between these symmetry functions.
Because of the absence of this non-linear interaction of symmetry functions, there is no clear
separation between all phases in the original feature vector (see figure 5.12d). Hidden layer 1
132
considers this non-linear interaction between the symmetry function, as result of which it is able
to understand the representation of 2H and 1T’ crystal structure (figure 5.12e). However, it is still
not able to distinguish between defects and 1T’ crystal structure. Addition of more layers alleviate
this problem, and the third layer of the NN model is able to distinguish all the six phases (see figure
5.12-f).
5.4.4 Neural network training and feature learning
The training accuracy of the model after 200 epochs is 93 %, where each atom is
represented by a 436-dimension vector. However, in this model, a large number feature vector
which are radial and angular symmetry function are used to represent the local environment of the
atom. Since, the hyper-parameters of the symmetry function are chosen randomly, hence, many
attributes of the feature vector will be highly correlated or irrelevant. Removing these irrelevant
features will not only make the computation faster but also it will make the model more compact
and increase its interpretability.
There are many ways to extract the relevant features directly from the model. For example,
in linear regression model we can look at the weights of feature vector in the trained model so as
to find out their importance. Here, features with higher weight are important as compared to the
ones with low weight. Another approach is based on the sensitivity analysis in which a particular
feature is removed from the feature vector and the model is retrained, and we do that for all the
features. Here, the feature that shows the least change in training accuracy is the irrelevant feature.
Hence to remove all the irrelevant features, this process is repeated until the training accuracy
decreases by a significant amount. However, it is very time consuming and not feasible when the
feature dimension is very large.
133

Figure 5.13: (a) Shows 10 features whose gradient change is maximum during training using full
feature vector set. (b) 10 features whose gradient change is minimum during training. (c) and (d)
shows training and validation accuracy of models, which are training with 50 features with
maximum gradient change or minimum gradient change in the original feature space.
For multi-layer neural network, there is another approach to find the relevant features,
which is based on the gradient of predicted output with respect to the input feature. In case of linear
model, this gradient is constant, which is not the case for multi-layer neural network. This gradient
can be computed for each feature progressively during the training process and the components of
the features whose rate of change in gradient is maximum are the relevant ones. This happens
because the gradient of the relevant features with respect to the output changes more rapidly as
comparted to the irrelevant ones for tiny changes in the output or input values. Since in
classification task we have more than one output class, and the gradient can be positive or negative;
134
therefore, we take the summation of the absolute value of the gradient for all outputs of the entire
training data set as proposed by Ruck et al in 1990:
𝑆 𝑖 =∑ ∑
𝜕 𝑓 𝑗 (𝑥 𝑙 )
𝜕 𝑥 𝑖 𝑔 𝑗 =1
𝑁 𝑙 =1
Eq 5. 11
Here, N is the total number of training examples, g is the number of output classes, xi is
the specific feature whose gradient we are computing, xl is the l
th
training example and fj(x
l
) is the
output of the NN mode.
Figure 5.13a-b shows the gradient of the 10 most features whose gradient either changed
the most (figure 5.13a) or the least (figure 5.13b). We haven trained another model using only a
50-dimension feature, which consists of only the most significant features with maximum gradient
change or the least significant ones. The training accuracy of these two models is shown in figure
5.13c-d, where we can see that model constructed using only the 50 most significant features have
an accuracy of 91% which is just 2% less then model accuracy created using entire feature vector.
On the other hand, model constructed using less significant features have only a training accuracy
of 58%. Further analysis of the symmetry function present in the relevant and irrelevant feature
set shows that relevant feature set consist of mostly radial symmetry functions with few angular
features. Figure 5.14a-b shows the couple of symmetry functions observed in the relevant and
irrelevant feature set. For example, relevant feature set (RFS) contains angular feature for Mo-Se-
Mo bond with bond angle peak at 0
0
, which make sense because three of the projected bond angle
between the top and bottom layer Se atom with Mo is close to zero in 2H phase. Similarly, bond
length of Mo-Se is around 2.8 Å whereas in defects Mo-Se bond is more than 3Å distance, hence,
RFS captures the radial symmetry function which are centered around these peaks.
135

Figure 5.14: 10 most significant radial and angular symmetry functions learned by the neural
network.
Following analysis shows, a careful feature analysis will not only help us to remove
irrelevant features from the model but can also be used to underlying structure of the atomic
configuration the NN model has looked at so as to identify all the phases.
5.4.5 Hyper-parameter tuning and data visualization
Figure 5.15a visualization of the fracture data set where classes are predicted by the model.
We can observe that, it classifies the 1T and 2H crystal structure accurately; however, it mislabels
some of the defects as 1T and vice versa. Remember in the t-SNE visualization of the hidden layer
(figure 5.12f), clusters of 2H and 1T crystal structures are completely separated whereas the
boundary between the 1T phase and defects is not clear. This shows that for some atoms, their
local neighbor environment is similar in 1T phase and defects, and this causes this
misclassification. For example, for some defects Mo-Se-Mo bond angle (BA) does not deviated
much from the mean BA value of crystals as here Mo-Se bond length is slightly above bond length
of crystals, which is 2.9 Å. Since, in this model a cutoff distance of 8Å is used for symmetry
function calculation, and hence damping factor will be very small for these defects. This creates
noise or ambiguity in the angular symmetry function value as it will exactly same in the 1T phase
136
and the defects. Figure 5.15b shows the visualization of the same fracture data but here the cutoff
distance for angular symmetry function is reduced to 3Å. This tuning of the cutoff distance
increased the model accuracy up to 95% as we can observe that misclassification between 1T phase
and defects has significantly reduced.

Figure 5.15: Visualization of validation set of a MoWSe2 fracture data for modes tuned with two
different cutoffs for angular symmetry function calculation (a) rcut= 8 Å and (b) rcut=3 Å.
5.4.6 Summary
In summary, we observe neural network-based analysis code can automatically learn the
relevant information from the simulation data. A careful tuning various model hyperparameters
like the number of hidden layers, dimension and the of the radial and angular symmetry function,
cutoff distance for the symmetry function calculation, increases the model performance up to 95%
137
with very small model complexity. The analysis of the learned features of the network gives us
further information about the data like the similarity between different phases
5.5 Conclusion
In conclusion, MD simulations of mode 1 fracture in MoWSe2 monolayer reveal that crack
propagation in the alloyed system is accompanied by a structural phase transformation from the
2H to the 1T phase in the process zone around the crack. Crack propagation mechanism in these
systems is highly directional dependent, where (1,1) loading shows crack branching and crack
closure, whereas in (1,0) loading direction shows cleavage fracture. Toughing mechanism during
crack propagation in these monolayers includes structural phase transformation from the 2H to the
1T phase, crack branching, crack closure, secondary grain nucleation. Also, these alloyed system
shows higher toughness as comparted to pure systems due to lower crack speed and process zone
volume
An automated machine learning algorithm is developed to analyze this strain induced phase
transition in MoWSe2. The algorithm takes as only the radial and angular description of atoms as
input and classify them into different phases which are 2H, 1T’ and defects. The neural network
work even learns the relevant features directly during the training process which the complexity
of the model significantly.

138
Conclusion
From the beginning the civilization, exploration and design of material has given us different
types of materials, which are polymers, ceramics, metals, etc. Some of these materials are just
simply discovered by us while there is another class of materials that are designed by humans by
altering the underlying chemistry and the atomic arrangement inside the material. However, the
strength and application of each materials is limited because of their inherent flaws. For example,
ceramics have high strength, hardness but are very brittle in nature. We can learn a lot from nature
here as observe that biological material designs themselves in a way such that they can service any
harsh environment. For example, the shell of sea creatures is nothing but a ceramic material, yet
they have very high fracture toughness. In nature, all biological material exhibits two unique
properties that are hierarchy in their microstructure and self-healing phenomena. Structural
hierarchy makes biological material deformable even when its constituent building block are
ceramic and also give them high strength and toughness. Their self-healing that further increases
their reliability, damage tolerance and service life. Taking a cue from nature, material scientists
are designing materials who can these natural phenomena. These are classified into two classes
that are micro-architectured material and self-healing material.
Micro-architectured material uses structural hierarchy to design materials. Using hierarchical
structure, we can create materials with very low density and high strength or make a ceramic
behave like a ductile material. In general, their mechanical property depends on the composition
of the base material, cell shape and topology and the struct length and thickness. However, the
effect of these design choices on the deformation mechanism of the hierarchical structure is not
well understood. Such understanding will help us create new class of materials. Hence, in this
work, we have studied the deformation mechanism of solid and hollow Ni kagomé lattices. Our
139
molecular dynamics revels that that hollow kagomé unit cell (HKUC) have higher strength and
toughness as compared to solid kagomé unit cell (SKUC). We observe that HKUC have higher
yield point (3.9 % strain) compared to SKUC (2.5% strain). Further, HKUCs shows two stage
deformation process where up to 11% strain, the deformation is highly localized near the nodal
region and for strain above 11%, all the struts and node starts deforming. Due to this localized
deformation in HKUC up to 11% strain, the maximum load bearing capacity of the system does
not drop much, and system shows less bending of its struts.
In the second project, we looked materials design using self-healing phenomena. In particular,
we looked at the self-healing phenomena of ceramics. Ceramics are lightweight, have high strength
and stiffness, and can withstand high temperatures. However, technological applications of
ceramics are limited by brittleness and sensitivity to flaws and cracks. In this work, we have
examined the atomistic mechanisms of crack healing in Al2O3 by oxidized n-SiC. Our reactive
molecular dynamics (RMD) simulations of Al2O3 containing SiC/SiO2 nanoparticles show that
silica is an excellent healing agent of cracks in alumina at high temperature. When the applied
strain is increased, bonds break and nanovoids are formed in the alumina matrix between the crack
front and nanoparticles. These nanovoids are the conduits through which silica diffuses into the
crack to stop its growth. The self-diffusion coefficient of silica is liquid-like (> 10
-5
cm
2
/s), and it
increases rapidly with an increase in the applied strain. As a result, the higher the applied strain,
the more rapid the crack healing. RMD simulations also show that alumina cavities containing
SiC/SiO2 nanoparticle are faceted, and the low energy prismatic (A) {2
̅
110} and prismatic (M)
{1
̅
010} planes are the prominent facets. These results are consistent with TEM studies of
nanocavities in -Al2O3. We observe grain formation on facets and grain growth with an increase
in the strain. Damage in the form of nanovoids and nanocracks is observed in grain boundaries,
140
but rapid diffusion of silica from nanoparticles again heals the damage. The weak link in the
strained Al2O3 matrix is the [011
̅
1] slip direction in the prism plane, along which nanovoids can
nucleate and coalesce to form cracks and cause fracture in the [011
̅
1] direction. Nanovoid growth
can be mitigated by a judicious choice of the size and spatial distribution of SiC/SiO2 nanoparticles
along the slip direction.
In the third project, we looked at a phase transition based toughening mechanism in
materials. In this toughening mechanism, volume change associated with the transformation causes
crack closure. Here, we have looked at the mechanism of crack propagation and its directional
dependence on the loading direction in MoWSe2 monolayer. Our MD simulations of mode 1
fracture in MoWSe2 monolayer reveal that crack growth in these alloyed systems is accompanied
by a structural phase transformation from the 2H to the 1T’ phase in the process zone around the
crack. Crack propagation mechanism in these systems is highly directional dependent, where (1,1)
loading shows crack branching and crack closure, whereas in (1,0) loading it shows cleavage
fracture. Toughing mechanism during crack propagation includes structural phase transformation
from the 2H to the 1T phase, crack branching, crack closure, secondary grain nucleation. Also,
these alloyed system shows higher toughness as comparted to pure systems due to lower crack
speed and process zone volume.
We have further developed an automated machine learning algorithm to analyze this strain
induced phase transition in MoWSe2. The algorithm takes as only the radial and angular
description of atoms as input and classify them into different phases which are 2H, 1T’ and defects.
The neural network work learns the relevant features directly during the training process which
the complexity of the model significantly.
141
References
1 Sottos, N., White, S. & Bond, I. Introduction: self-healing polymers and composites.
Journal of The Royal Society Interface 4, 347-348, (2007).
2 Lakes, R. Materials with structural hierarchy. Nature 361, 511, (1993).
3 Fleck, N. A., Deshpande, V. S. & Ashby, M. F. Micro-architectured materials: past,
present and future. Proceedings of the Royal Society A: Mathematical, Physical and
Engineering Science 466, 2495 (2010).
4 Schaedler, T. A. et al. Ultralight Metallic Microlattices. Science 334, 962-965, (2011).
5 Bauer, J., Hengsbach, S., Tesari, I., Schwaiger, R. & Kraft, O. High-strength cellular
ceramic composites with 3D microarchitecture. Proceedings of the National Academy of
Sciences 111, 2453 (2014).
6 Meza, L. R., Das, S. & Greer, J. R. Strong, lightweight, and recoverable three-
dimensional ceramic nanolattices. Science 345, 1322 (2014).
7 Montemayor, L. C. & Greer, J. R. Mechanical Response of Hollow Metallic Nanolattices:
Combining Structural and Material Size Effects. Journal of Applied Mechanics 82,
071012-071010, (2015).
8 Maggi, A., Li, H. & Greer, J. R. Three-dimensional nano-architected scaffolds with
tunable stiffness for efficient bone tissue growth. Acta Biomaterialia 63, 294-305, (2017).
9 Dooley, C. & Taylor, D. Self‐healing materials: what can nature teach us? Fatigue &
Fracture of Engineering Materials & Structures 40, 655-669, (2017).
10 Bekas, D. G., Tsirka, K., Baltzis, D. & Paipetis, A. S. Self-healing materials: A review of
advances in materials, evaluation, characterization and monitoring techniques.
Composites Part B: Engineering 87, 92-119, (2016).
11 Hayes, S. A., Zhang, W., Branthwaite, M. & Jones, F. R. Self-healing of damage in fibre-
reinforced polymer-matrix composites. Journal of the Royal Society Interface 4, 381-387,
(2007).
12 Williams, K. A., Boydston, A. J. & Bielawski, C. W. Towards electrically conductive,
self-healing materials. Journal of the Royal Society Interface 4, 359-362, (2007).
13 Verberg, R., Dale, A. T., Kumar, P., Alexeev, A. & Balazs, A. C. Healing substrates with
mobile, particle-filled microcapsules: designing a ‘repair and go’ system. Journal of The
Royal Society Interface 4, 349-357, (2007).
14 Kalista, S. J. & Ward, T. C. Thermal characteristics of the self-healing response in
poly(ethylene-co-methacrylic acid) copolymers. Journal of the Royal Society Interface 4,
405-411, (2007).
15 Madhan, M. & Prabhakaran, G. Self-healing Ability of Structural Ceramics – A Review.
Trends in Intelligent Robotics, Automation, and Manufacturing, 466-474, (2012).
16 Takahashi, K., Ando, K., Murase, H., Nakayama, S. & Saito, S. Threshold Stress for
Crack-Healing of Si3N4/SiC and Resultant Cyclic Fatigue Strength at the Healing
Temperature. Journal of the American Ceramic Society 88, 645-651, (2005).
142
17 Ando, K., Chu, M. C., Matsushita, S. & Sato, S. Effect of crack-healing and proof-testing
procedures on fatigue strength and reliability of Si3N4/SiC composites. Journal of the
European Ceramic Society 23, 977-984, (2003).
18 Thompson, A. M., Chan, H. M., Harmer, M. P. & Cook, R. E. Crack Healing and Stress
Relaxation in Al2O3 SiC “Nanocomposites”. Journal of the American Ceramic Society
78, 567-571, (1995).
19 Zhao, J. et al. Mechanical Behavior of Alumina–Silicon Carbide „Nanocomposites”.
Journal of the American Ceramic Society 76, 503-510, (1993).
20 Li, S., Song, G., Kwakernaak, K., van der Zwaag, S. & Sloof, W. G. Multiple crack
healing of a Ti2AlC ceramic. Journal of the European Ceramic Society 32, 1813-1820,
(2012).
21 Song, G. M. et al. Oxidation-induced crack healing in Ti3AlC2 ceramics. Scripta
Materialia 58, 13-16, (2008).
22 Novoselov, K. S. et al. Two-dimensional atomic crystals. Proceedings of the National
Academy of Sciences of the United States of America 102, 10451 (2005).
23 Feng, Y. et al. Synthesis of Large-Area Highly Crystalline Monolayer Molybdenum
Disulfide with Tunable Grain Size in a H2 Atmosphere. ACS Applied Materials &
Interfaces 7, 22587-22593, (2015).
24 Butler, S. Z. et al. Progress, Challenges, and Opportunities in Two-Dimensional
Materials Beyond Graphene. ACS Nano 7, 2898-2926, (2013).
25 Wang, Q. H., Kalantar-Zadeh, K., Kis, A., Coleman, J. N. & Strano, M. S. Electronics
and optoelectronics of two-dimensional transition metal dichalcogenides. Nature
Nanotechnology 7, 699, (2012).
26 Radisavljevic, B., Radenovic, A., Brivio, J., Giacometti, V. & Kis, A. Single-layer MoS2
transistors. Nature Nanotechnology 6, 147, 2011.
27 Kang, K. et al. High-mobility three-atom-thick semiconducting films with wafer-scale
homogeneity. Nature 520, 656, (2015).
28 Meza, L. R. et al. Reexamining the mechanical property space of three-dimensional
lattice architectures. Acta Materialia 140, 424-432, (2017).
29 Ashby, M. F. Cellular solids - scaling of properties. Cellular Ceramics, 3-17, (2005).
30 Deshpande, V. S., Ashby, M. F. & Fleck, N. A. Foam topology: bending versus
stretching dominated architectures. Acta Materialia 49, 1035-1040, (2001).
31 Wallach, J. C. & Gibson, L. J. Mechanical behavior of a three-dimensional truss material.
International Journal of Solids and Structures 38, 7181-7196, (2001).
32 Zhang, Y. H., Qiu, X. M. & Fang, D. N. Mechanical Properties of two novel planar
lattice structures. International Journal of Solids and Structures 45, 3751-3768, (2008)
33 Montemayor, L. C., Wong, W. H., Zhang, Y. W. & Greer, J. R. Insensitivity to Flaws
Leads to Damage Tolerance in Brittle Architected Meta-Materials. Scientific Reports 6,
20570, (2016).
143
34 Montemayor Lauren, C., Meza Lucas, R. & Greer Julia, R. Design and Fabrication of
Hollow Rigid Nanolattices via Two‐Photon Lithography. Advanced Engineering
Materials 16, 184-189, (2013).
35 Jacobsen, A. J., Barvosa-Carter, W. & Nutt, S. Compression behavior of micro-scale
truss structures formed from self-propagating polymer waveguides. Acta Materialia 55,
6724-6733, (2007).
36 Jiang, Y. & Wang, Q. Highly-stretchable 3D-architected Mechanical Metamaterials.
Scientific Reports 6, 34147, (2016).
37 Daw, M. S., Foiles, S. M. & Baskes, M. I. The embedded-atom method: a review of
theory and applications. Materials Science Reports 9, 251-310, (1993).
38 Uchic, M. D., Shade, P. A. & Dimiduk, D. M. Micro-compression testing of fcc metals:
A selected overview of experiments and simulations. Jom 61, 36-41, (2009).
39 Mishin, Y., Farkas, D., Mehl, M. J. & Papaconstantopoulos, D. A. Interatomic potentials
for monoatomic metals from experimental data and ab initio calculations. Physical
Review B 59, 3393-3407, (1999).
40 Tsuzuki, H., Branicio, P. S. & Rino, J. P. Structural characterization of deformed crystals
by analysis of common atomic neighborhood. Computer Physics Communications 177,
518-523, (2007).
41 Stukowski, A. & Albe, K. Extracting dislocations and non-dislocation crystal defects
from atomistic simulation data. Modelling and Simulation in Materials Science and
Engineering 18, 085001, (2010)
42 Wen, Y. H., Zhu, Z. Z., Shao, G. F. & Zhu, R. Z. The uniaxial tensile deformation of Ni
nanowire: atomic-scale computer simulations. Physica E-Low-Dimensional Systems &
Nanostructures 27, 113-120, (2005).
43 Setoodeh, A. R., Attariani, H. & Khosrownejad, M. Nickel nanowires under uniaxial
loads: A molecular dynamics simulation study. Computational Materials Science 44,
378-384, (2008).
44 Peng, C. et al. Strain rate dependent mechanical properties in single crystal nickel
nanowires. Applied Physics Letters 102, 083102, (2013).
45 Setoodeh, A. R. & Attariani, H. Nanoscale simulations of Bauschinger effects on a nickel
nanowire. Materials Letters 62, 4266-4268, (2008).
46 Branicio, P. S. & Rino, J. P. Large deformation and amorphization of Ni nanowires under
uniaxial strain: A molecular dynamics study. Physical Review B 62, 16950-16955,
(2000).
47 Lian, J. et al. Catastrophic vs Gradual Collapse of Thin-Walled Nanocrystalline Ni
Hollow Cylinders As Building Blocks of Microlattice Structures. Nano Letters 11, 4118-
4125, (2011).
48 Rabkin, E., Nam, H. S. & Srolovitz, D. J. Atomistic simulation of the deformation of gold
nanopillars. Acta Materialia 55, 2085-2099, (2007).
144
49 Greer, J. R. & Nix, W. D. Nanoscale gold pillars strengthened through dislocation
starvation. Physical Review B 73, 245410, (2006).
50 Park, H. S., Gall, K. & Zimmerman, J. A. Deformation of FCC nanowires by twinning
and slip. Journal of the Mechanics and Physics of Solids 54, 1862-1881, (2006).
51 Afanasyev, K. A. & Sansoz, F. Strengthening in gold nanopillars with nanoscale twins.
Nano Letters 7, 2056-2062, (2007).
52 Dongare, A. M., LaMattina, B. & Rajendran, A. M. Strengthening Behavior and
Tension–Compression Strength–Asymmetry in Nanocrystalline Metal–Ceramic
Composites. Journal of Engineering Materials and Technology 134, 041003-041003-
041008, (2012).
53 Neel, C. & Thadhani, N. Shock and release wave speed of an alumina epoxy composite.
Journal of Applied Physics 106, 046105, (2009).
54 Sinnott, S. B. & Dickey, E. C. Ceramic/metal interface structures and their relationship to
atomic- and meso-scale properties. Materials Science and Engineering: R: Reports 43, 1-
59, (2003).
55 Mauldin, T. C. & Kessler, M. R. Self-healing polymers and composites. International
Materials Reviews 55, 317-346, (2010).
56 Lumley, R. N., Morton, A. J. & Polmear, I. J. Enhanced creep performance in an Al–Cu–
Mg–Ag alloy through underageing. Acta Materialia 50, 3597-3608, (2002).
57 Zhang, S. et al. Self Healing of Creep Damage by Gold Precipitation in Iron Alloys.
Advanced Engineering Materials 17, 598-603, doi:10.1002/adem.201400511 (2015).
58 Swallow, J. G. et al. Chemomechanics of ionically conductive ceramics for electrical
energy conversion and storage. Journal of Electroceramics 32, 3-27, (2014).
59 Chroneos, A., Yildiz, B., Tarancon, A., Parfitt, D. & Kilner, J. A. Oxygen diffusion in
solid oxide fuel cell cathode and electrolyte materials: mechanistic insights from
atomistic simulations. Energy & Environmental Science 4, 2774-2789, (2011).
60 Coillot, D., Méar, F. O., Podor, R. & Montagne, L. Autonomic Self-Repairing Glassy
Materials. Advanced Functional Materials 20, 4371-4374, (2010).
61 Coillot, D., Méar, F. O., Podor, R. & Montagne, L. Influence of the Active Particles on
the Self-Healing Efficiency in Glassy Matrix. Advanced Engineering Materials 13, 426-
435, (2011).
62 Ando, K., Kim, B. S., Chu, M. C., Saito, S. & Takahashi, K. Crack-healing and
mechanical behaviour of Al2O3/SiC composites at elevated temperature. Fatigue &
Fracture of Engineering Materials & Structures 27, 533-541, (2004).
63 Nakao, W., Ono, M., Lee, S.-K., Takahashi, K. & Ando, K. Critical crack-healing
condition for SiC whisker reinforced alumina under stress. Journal of the European
Ceramic Society 25, 3649-3655, (2005).
64 Osada, T., Nakao, W., Takahashi, K., Ando, K. & Saito, S. Strength recovery behavior of
machined Al2O3/SiC nano-composite ceramics by crack-healing. Journal of the
European Ceramic Society 27, 3261-3267, (2007).
145
65 Osada, T., Nakao, W., Takahashi, K. & Ando, K. Kinetics of Self-Crack-Healing of
Alumina/Silicon Carbide Composite Including Oxygen Partial Pressure Effect. Journal of
the American Ceramic Society 92, 864-869, (2009).
66 Nihara, K. New Design Concept of Structural Ceramics - Ceramic Nanocomposites.
Nippon Seram Kyo Gak 99, 974-982 (1991).
67 Vashishta, P., Kalia, R. K., Nakano, A. & Rino, J. P. Interaction potentials for alumina
and molecular dynamics simulations of amorphous and liquid alumina. Journal of
Applied Physics 103, 083504, (2008).
68 Vashishta, P., Kalia, R. K., Nakano, A. & Rino, J. P. Interaction potential for silicon
carbide: A molecular dynamics study of elastic constants and vibrational density of states
for crystalline and amorphous silicon carbide. Journal of Applied Physics 101, 103515,
(2007).
69 Vashishta, P., Kalia, R. K., Rino, J. P. & Ebbsjo, I. Interaction Potential for Sio2: a
Molecular-Dynamics Study of Structural Correlations. Physical Review B 41, 12197-
12209, (1990).
70 Choi, J.-H. et al. Equilibrium Shape of Internal Cavities in Sapphire. Journal of the
American Ceramic Society 80, 62-68, (1997).
71 Mullins, W. W. & Rohrer, G. S. Nucleation Barrier for Volume-Conserving Shape
Changes of Faceted Crystals. Journal of the American Ceramic Society 83, 214-216,
(2000).
72 Kaplan, W. D., Levin, I. & Brandon, D. G. Significance of faceting on SiC nanoparticles
in alumina. Mater Sci Forum 207-209, 733-736, (1996).
73 Gupta, S., Zhang, Q., Emrick, T., Balazs, A. C. & Russell, T. P. Entropy-driven
segregation of nanoparticles to cracks in multilayered composite polymer structures.
Nature Materials 5, 229, (2006).
74 Upmanyu, M., Srolovitz, D. J., Lobkovsky, A. E., Warren, J. A. & Carter, W. C.
Simultaneous grain boundary migration and grain rotation. Acta Materialia 54, 1707-
1719, (2006).
75 Wang, Y. B., Li, B. Q., Sui, M. L. & Mao, S. X. Deformation-induced grain rotation and
growth in nanocrystalline Ni. Applied Physics Letters 92, 011903, (2008).
76 Thompson, A. P., Plimpton, S. J. & Mattson, W. General formulation of pressure and
stress tensor for arbitrary many-body interaction potentials under periodic boundary
conditions. The Journal of Chemical Physics 131, 154107, (2009).
77 Falk, M. L. Molecular-dynamics study of ductile and brittle fracture in model
noncrystalline solids. Physical Review B 60, 7062-7070, (1999).
78 Zhang, C., Kalia, R. K., Nakano, A., Vashishta, P. & Branicio, P. S. Deformation
mechanisms and damage in α-alumina under hypervelocity impact loading. Journal of
Applied Physics 103, 083508, (2008).
79 Heuer, A. H., Lagerlöf, K. P. D. & Castaing, J. Slip and twinning dislocations in sapphire
(α-Al2O3). Philosophical Magazine A 78, 747-763, (1998).
146
80 Lagerlöf, K. P. D., Heuer, A. H., Castaing, J., Rivière, J. P. & Mitchell, T. E. Slip and
Twinning in Sapphire (α-Al2O3). Journal of the American Ceramic Society 77, 385-397,
(1994).
81 Gupta, A., Sakthivel, T. & Seal, S. Recent development in 2D materials beyond
graphene. Progress in Materials Science 73, 44-126, (2015).
82 Venkata Subbaiah, Y., Saji, K. & Tiwari, A. Atomically thin MoS2: a versatile
nongraphene 2D material. Advanced Functional Materials 26, 2046-2069 (2016).
83 Zhang, Y., Tan, Y.-W., Stormer, H. L. & Kim, P. Experimental observation of quantum
Hall effect and Berry's phase in graphene. arXiv preprint cond-mat/0509355 (2005).
84 Liu, K. et al. Elastic Properties of Chemical-Vapor-Deposited Monolayer MoS2, WS2,
and Their Bilayer Heterostructures. Nano Letters 14, 5097-5103, (2014).
85 Island, J. O. et al. Precise and reversible band gap tuning in single-layer MoSe2 by
uniaxial strain. Nanoscale 8, 2589-2593, (2016).
86 Lee, Y. H. et al. Synthesis of Large‐Area MoS2 Atomic Layers with Chemical Vapor
Deposition. Advanced Materials 24, 2320-2325, (2012).
87 Ambrosi, A., Sofer, Z. & Pumera, M. 2H →1T phase transition and hydrogen evolution
activity of MoS2, MoSe2, WS2 and WSe2 strongly depends on the MX2 composition.
Chemical Communications 51, 8450-8453, (2015).
88 Li, H., Wu, J., Yin, Z. & Zhang, H. Preparation and Applications of Mechanically
Exfoliated Single-Layer and Multilayer MoS2 and WSe2 Nanosheets. Accounts of
Chemical Research 47, 1067-1075, (2014).
89 Apte, A. et al. Structural Phase Transformation in Strained Monolayer MoWSe2 Alloy.
ACS Nano 12, 3468-3476, (2018).
90 Lee, C. et al. Anomalous Lattice Vibrations of Single- and Few-Layer MoS2. ACS Nano
4, 2695-2700, (2010).
91 Krishnamoorthy, A. et al. Semiconductor-metal structural phase transformation in MoTe2
monolayers by electronic excitation. Nanoscale 10, 2742-2747, (2018).
92 Stillinger, F. H. & Weber, T. A. Computer simulation of local order in condensed phases
of silicon. Physical Review B 31, 5262-5271, (1985).
93 Roy, A. et al. Structural and Electrical Properties of MoTe2 and MoSe2 Grown by
Molecular Beam Epitaxy. ACS Applied Materials & Interfaces 8, 7396-7402, (2016).
94 Zhang, R., Koutsos, V. & Cheung, R. Elastic properties of suspended multilayer WSe2.
Applied Physics Letters 108, 042104, (2016).
95 Yang, Y. et al. Brittle Fracture of 2D MoSe2. Advanced Materials 29, 1604201, (2016).
96 Rajan, K. Materials informatics. Materials Today 8, 38-45 (2005).
97 LeSar, R. Materials informatics: An emerging technology for materials development.
Statistical Analysis and Data Mining: The ASA Data Science Journal 1, 372-374 (2009).
98 Mueller, T., Kusne, A. G. & Ramprasad, R. Machine learning in materials science:
Recent progress and emerging applications. Rev. Comput. Chem (2015).
147
99 Lee, J., Seko, A., Shitara, K., Nakayama, K. & Tanaka, I. Prediction model of band gap
for inorganic compounds by combination of density functional theory calculations and
machine learning techniques. Physical Review B 93, 115104 (2016).
100 Gu, T., Lu, W., Bao, X. & Chen, N. Using support vector regression for the prediction of
the band gap and melting point of binary and ternary compound semiconductors. Solid
state sciences 8, 129-136 (2006).
101 Kim, C., Pilania, G. & Ramprasad, R. From organized high-throughput data to
phenomenological theory using machine learning: the example of dielectric breakdown.
Chemistry of Materials 28, 1304-1311 (2016).
102 Kim, C., Pilania, G. & Ramprasad, R. Machine learning assisted predictions of intrinsic
dielectric breakdown strength of ABX3 perovskites. The Journal of Physical Chemistry C
120, 14575-14580 (2016).
103 Zhaochun, Z., Ruiwu, P. & Nianyi, C. Artificial neural network prediction of the band
gap and melting point of binary and ternary compound semiconductors. Materials
Science and Engineering: B 54, 149-152, (1998).
104 Huan, T. D., Mannodi-Kanakkithodi, A. & Ramprasad, R. Accelerated materials property
predictions and design using motif-based fingerprints. Physical Review B 92, 014106
(2015).
105 Pilania, G., Gubernatis, J. E. & Lookman, T. Multi-fidelity machine learning models for
accurate bandgap predictions of solids. Computational Materials Science 129, 156-163
(2017).
106 Mannodi-Kanakkithodi, A., Pilania, G., Huan, T. D., Lookman, T. & Ramprasad, R.
Machine learning strategy for accelerated design of polymer dielectrics. Scientific reports
6, 20952 (2016).

Abstract (if available)

Abstract Nanostructured materials have exceptional mechanical properties and have a wide variety of applications ranging from manufacturing, aerospace to healthcare. However, the service life of these materials is often limited due to the formation of cracks and defects that leads to their eventual failure, thus insights into atomistic mechanisms underlying these processes are valuable inputs for the design of next-generation nanostructured materials. In this research work, large-scale molecular dynamics simulation, combined with machine learning and data mining techniques is used to understand three basic problems related to nanostructured materials design with superior mechanical properties, which are: (1) ultra-light material design using hierarchical structures like kagomé lattices, (2) self-healing of cracks in alumina composite embedded with SiC/SiO₂ nanoparticles, and (3) toughing of MoWSe₂ monolayer during crack propagation by novel strain-induced structural phase transformation. ❧ Many natural and man-made materials exhibit hierarchical structures like for example bird wings, Eiffel tower, Garabit Viaduct. This hierarchy in structures can be used for materials design at the atomic level as it allows us to tune several materials properties like stiffness, compressive strength, fracture toughness by just changing the underlying structural arrangement of materials. In the first work, we have studied the deformation behavior of an ultralight architecture consisting of hollow Ni nanotubes and solid nanorods arranged as a 3-D kagomé structure. Our result shows that hollow kagomé structure have higher strength and stiffness in comparison of solid kagomé structures. The hollow kagomé structure shows two stage deformation mechanism where up to 11% strain, deformation is highly localized near the nodal region and at higher strain they show bending-dominated deformation. Whereas in solid kagomé structures, large defamations are observed in the entire structure even at lower strain. ❧ Designing material systems that can sense and repair damage is of great interest. Inspired by biological systems, material scientists are developing self-healing materials in which damage initiation triggers an automatic release and transport of healing agent to the damage site and repairs it. In the second work, we have investigated crack healing in a ceramic nanocomposite of alumina containing SiC/SiO₂ nanoparticles. Our result shows that the nanoparticles closest to the advancing crack in alumina distort to create nanochannels through which silica flows toward the crack and heals the damage. Under applied strain, alumina matrix also forms facets at interfaces. These facets act as nucleation sites for multiple secondary amorphous grains in alumina matrix. Nanovoids form along the grain boundaries but these nanovoids are again healed by molten silica released from the nanoparticles. The mechanisms of crack healing in alumina will help experimentalists to design novel self-healing ceramics for a broad class of energy technologies. ❧ Two-dimensional materials exhibit different mechanical property in comparison of their bulk counterpart owing to their monolayer atomic thickness. In fact, properties of these layered materials can be tuned by mechanical deformation. In the third work, we have examined the atomistic mechanisms of defect formation and crack propagation in a strained monolayer MoWSe₂ alloy using molecular dynamics. Our result shows that under strain, the crack meanders as it propagates through the system and branches into daughter cracks when it crosses the MoSe₂/WSe₂ interface and enters the WSe₂ patches. The most dramatic change occurs in the process zone around the crack tip, where large stress concentration and the resulting biaxial tensile strain trigger an irreversible local structural phase transformation from the ground state 2H crystal structure to the 1T crystal structure. This phenomenon of crack branching, crack closure and strain induced phase transformation in MoWSe₂ alloy increases its fracture toughness.

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

Neural network for molecular dynamics simulation and design of 2D materials

PDF

High throughput computational framework for synthesis and accelerated discovery of dielectric polymer materials using polarizable reactive molecular dynamics and graph neural networks

PDF

Controlling electronic properties of two-dimensional quantum materials: simulation at the nexus of the classical and quantum computing eras

PDF

Nanomaterials under extreme environments: a study of structural and dynamic properties using reactive molecular dynamics simulations

PDF

Molecular dynamics study of energetic materials

PDF

Quantum molecular dynamics and machine learning for non-equilibrium processes in quantum materials

PDF

Atomistic modeling of the mechanical properties of metallic glasses, nanoglasses, and their composites

PDF

Modeling of bicontinuous nanoporous materials and investigation of their mechanical and hydraulic properties

Asset Metadata

Creator Rajak, Pankaj (author)

Core Title Large-scale molecular dynamics simulations of nano-structured materials

School Viterbi School of Engineering

Degree Doctor of Philosophy

Degree Program Materials Science

Publication Date 10/16/2018

Defense Date 06/21/2018

Publisher University of Southern California (original), University of Southern California. Libraries (digital)

Tag ceramic nano-composite,grain-growth,hierarchical meta-materials,kagome lattices,layered material,micro-architectured material,neural network,OAI-PMH Harvest,reactive molecular dynamics,self-healing ceramics,strain-induced phase transition

Format application/pdf (imt)

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Vashishta, Priya (committee chair), Branicio, Paulo (committee member), Kalia, Rajiv (committee member), Nakano, Aiichiro (committee member)

Creator Email pankajr127@gmail.com,rajak@usc.edu

Permanent Link (DOI) https://doi.org/10.25549/usctheses-c89-80832

Unique identifier UC11671952

Identifier etd-RajakPanka-6854.pdf (filename),usctheses-c89-80832 (legacy record id)

Legacy Identifier etd-RajakPanka-6854.pdf

Dmrecord 80832

Document Type Dissertation

Format application/pdf (imt)

Rights Rajak, Pankaj

Type texts

Source University of Southern California (contributing entity), University of Southern California Dissertations and Theses (collection)

Access Conditions The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the a...

Repository Name University of Southern California Digital Library

Repository Location USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA