Fracture and Wrinkling in Graphene with Topological Defects by Teng Zhang B. Eng. Engineering Mechanics, Dalian University of Technology, 2007 M. Eng. Solid Mechanics, Dalian University of Technology, 2010 A dissertation submitted in partial fulfillment of the requirements for the Degree of Doctor of Philosophy in the School of Engineering at Brown University Providence, Rhode Island May 2015 © Copyright 2015 by Teng Zhang This dissertation by Teng Zhang is accepted in its present form by the School of Engineering as satisfying the dissertation requirement for the degree of Doctor of Philosophy. Date_________________ _________________________________ Prof. Huajian Gao, Advisor Recommended to the Graduate Council Date_________________ _________________________________ Prof. Kyung-Suk Kim, Reader Date_________________ _________________________________ Prof. Brian Sheldon, Reader Approved by the Graduate Council Date_________________ _________________________________ Prof. Peter M. Weber Dean of the Graduate School iii Vita Teng Zhang was born on September 02, 1985 in Kaifeng, Henan Province, China. He attended Dalian University of Technology in 2003, and received the B. Eng. in Engineering Mechanics in 2007. Upon graduation, he was enrolled by Dalian University of Technology with entrance examination waiver. He was rewarded the M. Eng in Solid Mechanics in 2010. Subsequently, he entered the Solid Mechanics program at Brown University. His research focused on modeling and simulations of mechanical properties of graphene with topological defects. iv Acknowledgments First, I would like to thank my advisor, Prof. Huajian Gao, for accepting me as his student and offering me his kind support as well as valuable guidance throughout all my years at Brown University. He taught me how to be a capable scholar through his great personality, deep physical insights, and intensive focus on research. Working with him during the past four and a half years has been a wonderful experience which will also benefit me in my future research career. Second, I would like to express my appreciation to the members of my thesis committee, Prof. Kyung-Suk Kim and Prof. Brian Sheldon, for taking time to read my thesis and provide me with many valuable suggestions. I am very impressed with Prof. Kim’s ability to do excellent work in both experiment and theory. I have benefited a lot from several collaborations with him as well as his course “Linear Elasticity”. I also feel very lucky to have participated in a research project led by Prof. Sheldon. Third, I want to express my gratitude to my wonderful friends, including many former and current group members, postdocs and graduate students at Brown. I will particularly thank Dr. Xiaoyan Li for his help and encouragement on my research. I would also like to thank the faculty of Brown Solid Mechanics group for their high- quality graduate courses. Finally, I would like to thank my family- my parents and my wife. Without their support and understanding, I would never have completed this thesis. v I dedicate this to my Dear Wife. vi Contents Signature Page iii Vita iv Acknowledgements v Table of Contents vii List of Tables x List of Figures xi 1 Introduction 1 1.1 Mechanical properties of graphene with topological defects ……..………….1 1.2 Methodology of molecular dynamics ..………………………………………..3 1.3 Outline of the thesis .....……………………………………......………………..6 2 Flaw Insensitive Fracture in Nanocrystalline Graphene 8 2.1 Introduction ……………………………………………………………………8 2.2 Concepts of flaw tolerance/insensitivity .…………………………………….9 2.3 MD Simulations for fracture of nanocrystalline graphene ……………..12 2.4 Summary ………..………………………………………………………….24 3 Crack Interacting with Dislocation and Grain Boundary 25 3.1 Introduction ……………………………………...……………………………25 vii 3.2 Dislocation core in graphene ……………………………...……………….….26 3.3 Crack interaction with a dislocation in graphene ……………...……………...28 3.4 Crack along a grain boundary in graphene …..…………………..………...…33 3.5 Summary ……………………………………………………………...……….40 4 Defects controlled wrinkling and topological design in graphene 43 4.1 Introduction ……………………………………………………………………43 4.2 Continuum model for graphene with topological defects …..………...….……45 4.3 Mathematical analogy between defects and inhomogeneous growth ..49 4.4 Numerical methods for simulation …………………………...………..………51 4.5 An isolated dislocation in graphene ……………………..…..………….…..…53 4.6 Defects controlled graphene structures ………..………..…..………..….……59 4.7 Summary ………..……………………………………………………………..63 5 A Phase Field Crystal Method for Multifunctional Curved Graphene Design with Topological Defects 65 5.1 Introduction ……………………………………………………………………65 5.2 Previous studies on curved graphene ……………………………………..…..66 5.3 Phase field crystal method ………………………………………..…………..71 5.4 FEM implementation of PFC model ……..…………………………………..76 5.5 Preliminary results on using PFC method in curved graphene design ……....78 5.6 Summary ………………………………………………………………………84 6 Conclusions and Perspectives 86 6.1 Concluding remarks ……..……………………………………………………86 viii 6.2 Outlook of future research ….…………………………………………………87 Bibliography …………………………………………………………………………...90 ix List of Tables Table 3.1. Simulated failure strains of dislocation interacting with a linear elastic crack versus a Dugdale cohesive crack from MD simulations and continuum models. ……..32 x List of Figures Figure 1.1. Graphene structures. (a) A picture of perfect graphene structure from Wikipedia (http://en.wikipedia.org/wiki/Graphene). (b) Polycrystalline graphene with grain boundaries [14]. (c) Controlled defects with ion irradiation [25]. ………………….2 Figure 1.2. Wrinkling in graphene with topological defects. (a) Atomistic simulations for wrinkling of grain boundary in graphene [33]. (b) Experiments and simulations show the wrinkling in graphene with a dislocation dipole [36]. …………………………………3 Figure 2.1. Theoretical modeling of flaw tolerance. (a) Schematic drawing of a surface crack in a mineral plate. (b) Comparison of the fracture strength of a cracked mineral platelet predicted from the Griffith criterion of fracture with the strength of a perfect crystal. (a)-(b) are both from reference [31]. ……………………………………………10 Figure 2.2. Theoretical modeling of flaw tolerance in a cracked strip. (a) A thin strip with a central crack. The strip width is 2H and the crack length is 2a. (b) A thin strip with edge cracks. The strip width is 2H and the crack length is a. (c) Theoretical analysis of flaw tolerance in the thin strip. (a)-(c) are all from reference [78]. ……….………………..12 Figure 2.3. Simulation sample of a nanocrystalline graphene nanostrip. (a) Contours of out-of-plane displacement after equilibration at 300 K. (b)-(d) Representative defect structures revealed by differences in potential energy, which include (b) pentagon - heptagon and pentagon-octagon pairs, (c) vacancy and (d) polygon. (Scale bar: 5 nm). .14 Figure 2.4. Deformation and fracture behaviors of a nanocrystalline graphene nanostrip with a center hole under longitudinal tension. (a) Typical tensile stress-strain curve. (b-d) xi Tensile stress contours (unit: GPa) at ε=8%, 9.1%, 9.25% and 9.35%, respectively. In (b) and (c), micro-cracks (marked by squares) initiate from topological defects on grain boundaries. In (d) and (e), micro-cracks coalesce to form a big crack, eventually leading to failure. Note that the fracture surface is nearly perpendicular to the loading direction. (Scale bar: 5 nm). ………………………………………………………………………..16 Figure 2.5. Strength variations of single- and nano-crystalline graphene nanostrips with the radius of an internal hole. (a) Strength of armchair and zigzag single-crystalline graphene.  EX is the experimental value of graphene strength reported in Ref. [4]. The reader is referred to Ref. [28] for strength calculated from a multiscale method coupling quantum, molecular and continuum statics and to Ref. [82] for strength from MD simulations of a square graphene sheet (55 nm2) with a central slit at room temperature, at strain rate of 5×108 /s. (b) Strength of nanocrystalline graphene. (c) Normalized strength of single- and nano-crystalline graphene nanostrips.  m denotes the strength from MD simulations,  th refers to the limiting strength of the flaw tolerant nanocrystalline graphene strip at the holed cross section. (d) Strengths of single- and nano-crystalline graphene strips with an ellipitic hole. ………………………………..19 Figure 2.6. Tensile stress contours (unit: GPa) near a hole of radius 3 nm in single- and nano-crystalline graphene nanostrips during deformation. (a) Armchair and (b) zigzag graphene strips at a strain of 4%. Nano-crystalline graphene strips (c) with and (d) without the hole at a strain of 8%. (Scale bar: 5nm). ……………………………………21 Figure 2.7. Flaw insensitive fracture in an edge-notched nc-graphene strip. (a) Typical tensile stress-strain curve of a pre-edge-notched nc-graphene strip. (b) Tensile stress contours (unit: GPa) of the strip at different strains. ……………………………………23 xii Figure 3.1 Dislocations dipole in graphene. (a), Molecular mechanics (MM) model for a dislocation dipole in a graphene square sheet (length L of 20 or 40 nm) with dipole distance d varying from 0.24 to 6 nm. (b), Potential energy for the dislocation dipole from MM simulations and classical edge dislocation model. (c), Atomic configuration of the dislocation dipole in graphene with d=26b. (d), Core of a heptagon-pentagon dislocation in graphene shown in (c). (Scale bar: 2 nm). …………………………..………………27 Figure 3.2 Crack interaction with a dislocation in graphene. (a), Stress-strain curves from MD simulations for a cracked graphene strip with/without a dislocation ahead of the crack tip. (b), Stress contours (yy) for a graphene strip with crack only (εyy=0.027). (c), A zoom-in snapshot around the crack tip in (b). (d), Stress contours for a crack interacting with a dislocation-1 (εyy=0.035). (e), A zoom-in snapshot around the crack tip in (d). (f), Stress contours for a crack interacting with a dislocation-2 (εyy=0.019). (g), A zoom-in snapshot around the crack tip in (f). (h), Stress contours for a crack interacting with a dislocation-2 (εyy=0.019). (i), A zoom-in snapshot around the crack tip in (h). The scale bars are 10 nm in (b), (d), (f) and (h). .………………………………………….…29 Figure 3.3 Cracking along a grain boundary in graphene. (a), Bi-crystal graphene strip with left crack. (b), A rightward crack along a zigzag grain boundary (θ=9.43). (c), A leftward crack along a zigzag grain boundary (θ=9.43). The scale bars are 10 nm in (a) and 1 nm in (b) and (c).. ……………………………………………………………...….34 Figure 3.4 MD simulations of a crack propagating along an armchair-tilted grain boundary with misorientation angle θ =9.43. (a), Stress-strain curves for the rightward and leftward cracks. (b) and (c), Crack propagation patterns. (d), (f) and (h), Snapshots for stress (σ22) distribution near the crack tip for the rightward crack under different xiii applied strains. (e), (g) and (i), Snapshots for stress (σ22) distribution near the crack tip for the leftward crack under different applied strains. The scale bars are 10 nm in (b-c) and 1 nm in (d-i). ……………………………………………………..……………………….35 Figure 3.5 Initial stress (σ22) distributions near the crack tip for a crack along an armchair-tilted grain boundary. (a) MD simulations of a rightward crack. (b), MD simulation of a leftward crack. (c), Continuum models for (a). (d), Continuum model for (b). The scale bar is 1 nm in (a-b).. ……………………………………..………………37 Figure 3.6 MD simulations of a crack propagating along an armchair-tilted grain boundary with misorientation angle θ =21.8. (a), Stress-strain curves for the rightward and leftward crack propagation. (b) and (c), Crack propagation patterns for the rightward and leftward cracks. (d) and (f), Snapshots of stress (σ22) distributions near the crack tip for the rightward crack under different applied strains. (e) and (g), Snapshots of stress distributions near the crack tip for the leftward crack under different applied strains. The blue lines in (b) and (f) represent the lowest energy surface. The scale bars are 10 nm in (b-c) and 1 nm in (d-g). ……………………………………………..…………………38 Figure 4.1. Atomic configuration of a dislocation dipole in a graphene sheet of dimension 20×20 nm2 (only part of the region around the dislocation dipole is shown here). (a) A perspective view of deformation around the dislocation dipole. (b) Top view of the dislocation dipole. (c) Bond structures around the dislocation core in 3D. (d) Bond structures around the dislocation core in 2D projection. The color represents the scale of the out-of-plane displacement in (a) and potential energy in (b-d), respectively. (Scale bar: 1 nm) ………………………………………………………………………………44 xiv Figure 4.2. Buckled shapes around disclinations in a circular graphene disk of radius 5 nm. (a) The triangular lattice model adopted in continuum simulations. (b) The calculated 3D buckled shapes around positive and negative disclinations in the circular graphene disk from the continuum model. (c) Corresponding 3D buckled shapes from atomistic simulations. The colors represent eigenstrain field with rc =0.5 Å and out-of-plane deformation in (b) and (c), respectively. ………………………………………………52 Figure 4.3. An isolated dislocation in a square graphene sheet. (a) Schematic illustration. (b) A magnified view of atomic structure of a dislocation (heptagon-pentagon dipole) in graphene. (c) The eigenstrain field of dislocation used in the continuum model. (Scale bar: 1 nm). ……………………………………………………………………………….53 Figure 4.4. Comparison of stress fields around an isolated dislocation calculated from continuum modelling and atomistic simulations. Contour plots of stress components 11,22, 12 from (a-c) continuum modelling based on the generalized von Karman equation; (d-f) 3D atomistic simulations; (g-i) 2D atomistic simulations (without out-of- plane deformation); (j-l) classical plane-stress solution of an edge dislocation. (Scale bar: 1 nm.) …………………………………………………...……………………………….55 Figure 4.5. Comparison of the out-of-plane displacement field calculated from continuum modelling based on the generalized von Karman equation and atomistic simulations. (a) Contour plots of the out-of-plane displacement field calculated from continuum modelling. (b) A zoom-in snapshot of deformation around the dislocation core in (a). (c) Corresponding contour plots of the out-of-plane displacement field calculated from atomistic simulations. (d) A zoom-in snapshot of deformation around the dislocation core in (c). The scale bars are 10 nm in (a) and (c), 1 nm in (b) and (d), respectively. ………56 xv Figure 4.6. Effect of an applied tensile strain on wrinkling around a dislocation. (a) Continuum solutions under uniaxial and biaxial strains. (b) Corresponding atomistic simulations. (Scale bar: 1 nm). ……………………………………………………….....57 Figure 4.7. A sinusoidal graphene ruga induced by a periodic array of disclinations from continuum modeling and atomistic simulations. (a) The continuum eigenstrain filed in the undeformed configuration with rc=0.5 Å and (b) the configuration of the sinusoidal graphene ruga predicted from continuum modelling. (c) Potential energy contour (top view) and (d) the configuration of the sinusoidal graphene ruga from atomistic simulations. (Scale bar: 1 nm) ……………………………………………………….....60 Figure 4.8. A catenoid graphene funnel from atomic simulations and continuum modeling. (a) Catnoid surfaces for c=0.5 and 0.7, from which zmax is determined by cosh(zmax/c)c=1. (b) The initial atomic configuration of 1018 carbon atoms on a catenoid surface (Kusumaatmaja and Wales [134]) with distributed dislocations. (c) The relaxed atomic structure of (b). The color represents the potential energy scale. The corresponding (d) initial and (e) final structures predicted from continuum modelling. The color represents the eigenstrain scale with rc=0.5 Å. (scale bar: 1 nm) …………....62 Figure 5.1. Curved carbon nanostructures. (a) Atomic structure of fullerene (C60) (http://en.wikipedia.org/wiki/Fullerene). Two-fold symmetry axis views of the cubic D (b) and P (c) of schwarzite, emphasizing differences in the layout of the seven-membered rings [151]. ……………………………………………………………….……………68 Figure 5.2. Three-dimensional nanoporous graphene. (a) Fabrication of nanoporous graphene with the use of nanoprous Ni. (b) The structure of nanoporous graphene with and without nanoporous Ni. Raman spectral indicates that the 3D nanoporous sample xvi consists of high-quality monolayered graphene. (c) Atomic structures of the nanoporous graphene. Topological defects can be found in the curved surface [155]. ……………69 Figure 5.3. Creating a curved carbon structure from a geometrical point of view. (a) Unit cell for a toroidal carbon nanotube (Ref. [141]). (b) Bead molecular structure for a toroidal tube. (c) Bead molecular structure for hierarchical fullerene. (b) and (c) from the online blog: http://thebeadedmolecules.blogspot.com/. (d) Constructing curved graphene sheet from the triangular lattice network via the voronoi scheme. (e) Complicated curved carbon structures based on the triangular lattice network. (d) and (e) are from Ref. [143]………………………………………………………..…………………………69 Figure 5.4. The phase field crystal method. (a) Differences and similarities between the phase field and phase field crystal methods (from Prof. Voorhees’s talk at NIST). (b) Phase diagram of the phase field crystal model (Ref. [158]). …………………..…..…73 Figure 5.5. Phase field crystal simulations on a curved surface. (a) PFC simulation of the classical Thomson problem [147]. (b) PFC simulation on the surface of the Stanford bunny [147]. (c) Defects pattern on the catenoid surface with PFC simulation [148]. (d) Defect dynamics on a sinusoidal surface using PFC [168].. …………………..………75 Figure 5.6. UFL code of FEM implementation for Poisson equation in FEniCS [170]. ..78 Figure 5.7. FEM simulations of typical PFC patterns. (a) A strip pattern for   0.05 . (b) a triangular pattern for   0.2 , and (c) a homogenous pattern for   0.6 . ….79 Figure 5.8. A flow chart for designing curved graphene sheet via PFC method. (a) The initial curved manifold. (b) A continuum triangular pattern on the curved surface generated by PFC method. (c) A discrete triangular lattice network according to the xvii triangular pattern. (d) The graphene atomic structure through voronoi construction from the triangular network. …………………………………………………………………..79 Figure 5.9. Uniaxial tensile simulation of a sinusoidal graphene sheet with molecular dynamics. (a) Stress-strain curves of the sinusoidal graphene with different unit cell size under the uniaxial tensile loading. (b)-(c) Top and front views of a snapshot of stress (σxx) contours under applied strain ε=0.15 for one unit cell with length L=4 nm. (d)-(e) Top and front views of a snapshot of stress contours under applied strain ε=0.15 for four basic unit cell with length L=16 nm. (The scale bars are 4 nm)………………………………81 Figure 5.10. Edge crack propagation in a nanostrip of sinusoidal and flat graphene. (a) Stress-strain curves for flat and sinusoidal graphene sheets. The dot-dashed line represents the result of a larger sinusoidal graphene strip with length of 160 nm and width of 80 nm. (b) A nanostrip of sinusoidal graphene with an edge crack, the color showing out-of plane deformation. (c) Crack tip and atomic orientation of zigzag, armchair and sinusoidal graphene. (d) A top view of stress (σyy) contours near the crack tip during the tensile deformation. (e) Contours of out-of-plane displacements from a perspective view near the crack tip during the tensile deformation. The scale bars are 10 nm and 1 nm in (b) and (d), respectively. ……….………………………………………………………82 xviii Chapter 1 Introduction 1.1. Mechanical properties of graphene with topological defects As a monolayer of sp2-bonded carbon atoms densely packed in a honeycomb crystal lattice [1] (as shown in Fig. 1.1a), graphene exhibits exceptional electrical, mechanical and optical properties [2-4], making it an ideal candidate for a vast range of applications, such as flexible electronics [5, 6], field effect transistors [1, 7, 8], nanomechanical resonators [9]. While it is still challenging to fabricate large pieces of pristine single- crystalline graphene, large-area polycrystalline graphene (polygraphene) (as shown in Fig. 1.1b) containing internal grain boundaries has been successfully synthesized via chemical vapor deposition (CVD) [5, 10]. Motivated by this development, considerable research has been dedicated to the effect of grain boundaries on the electrical [11-14], mechanical [14-22] and thermal properties [23, 24] of graphene. Recent studies (Fig. 1.1c) have also shown that the position and type of topological defects in graphene can be controlled with ion irradiation [25], which suggests a promising way to tune the physical properties of graphene via defect engineering [26]. a b c 1 Figure 1.1. Graphene structures. (a) A picture of perfect graphene structurefrom Wikipedia (http://en.wikipedia.org/wiki/Graphene). (b) Polycrystalline graphene with grain boundaries [14]. (c) Controlled defects with ion irradiation [25]. Fracture is among the most prominent concerns about the mechanical behaviors of polygraphene. It has been recently shown via atomistic and quantum simulations that the fracture strength of a bi-crystalline graphene is sensitive to the misorientation angle across the grain boundary [15, 16, 20, 22]. Experiments revealed that the strength of polygraphene is about 35 GPa [14, 17], which is larger than most of the engineering materials but substantially lower than the ideal strength (about 130 GPa) of single- crystalline graphene [4]. This is not so surprising since many forms of defects (e.g., vacancies, voids, micro-cracks and chemical impurities) can lead to reduction in the strength of graphene [15, 17, 27-30]. In spite of these impressive progress on the fracture behaviors of polycrystalline graphene, some fundamental issues remain to be fully elucidated. For example, could nanocrystalline graphene (nc-graphene), which contains a network of nanoscale internal grain boundaries, be made flaw tolerant [31] and if so, what would be the critical condition? It is already known that the grain boundary strength in graphene is generally weaker than that of perfect graphene. However, experiments and atomistic simulations [18] have shown that intra-granular fracture happens in polycrystalline graphene. An interesting open question is whether grain boundary can be tougher than the perfect graphene lattice? Due to its atomic scale thickness, the deformation energy in a free standing graphene sheet can be easily released through out-of-plane wrinkles, which have been predicted by atomistic simulations [32-35] and verified by carefully designed experiments[36, 37]. So 2 far, there is still the lack of a complete analysis of the wrinkling in graphene with defects from the continuum point of view. It has also been shown that the out-of-plane displacement in graphene can be used to tune the electrical [38-40] and mechanical properties of graphene, such as anisotropic friction [41] and tunable wettability through controllable crumpling [42] of graphene on substrate. These recent studies are calling for systematic investigations of topological defects induced wrinkling in graphene. b a Figure 1.2. Wrinkling in graphene with topological defects. (a) Atomistic simulations for wrinkling of grain boundary in graphene [33]. (b) Experiments and simulations show the wrinkling in graphene with a dislocation dipole [36]. 1.2 Methodology of molecular dynamics Most of simulations in this thesis will be carried out with molecular dynamics (MD). Therefore I will make a brief introduction of the MD simulation, interatomic potential and common ensembles used in the simulations. Molecular dynamics (MD) is a computational simulation of the physical movements of atoms and molecules [43], which have become a popular and standard tool to investigate the fundamental mechanisms of the mechanical properties of nanostructures. 3 For the classical molecular mechanics, the atoms and molecules interact with each other following Newton’s second law, miri  f i , (1.1) where mi is the mass of the ith atom, ri is the Cartsian coordinates, and fi is the interaction force applied on ith atom from other atoms, and “.” denotes the time derivative. The force is determined by U fi   , (1.2) ri where U is the potential energy. In classical molecular dynamics simulation, the potential energy is usually described with some empirical functions with tunable parameters to fit the experimental data and first principle calculations. Interatomic potential is the central part of the classical MD simulations and sometimes also the bottleneck of accurately predicting the physical properties of the structures due to the lack of accurate potential. There have been several popular atomic potentials for carbon-carbon atom interaction with the ability to describe bond break and forming, such as Tersoff [44], REBO (reactive empirical bond-order) [45], AIREBO (adaptive intermolecular REBO) [46], EDIP (environment-dependent interaction potential) [47], LCBOPII (long-range reactive bond-order potential) [48], ReaxFF (reactive force field) [49], COMB (Charge-Optimized Many-Body) [50], SED-REBO (screened environment-dependent reactive empirical bond-order ) [51], NN (Neural Network) potential [52]. Until now, AIREBO/REBO is still the most popular atomistic potential for simulation of carbon materials, as it can capture most of the physical properties of carbon 4 structures and is already implemented in the widely used open source MD package LAMMPS (Large-scale Atomic/Molecular Massively Parallel Simulator) [53]. Although NN, LCBOPII and SED-REBO potentials are shown to be more accurate, their applications are still limited due to its complication in implementation and unavailability in the public open source software. All the atomistic simulations in the thesis are performed with AIREBO potential. It should be noted that the smaller cutoff distance in the switching function of AIREBO potential should be taken in the range of 1.92~2.0 Angstrom to avoid a known non-physical post-hardening behavior [15, 16, 54], while remaining consistent with DFT calculations [20]. For the present study, this cutoff distance was set at 2.0 Å following a previous study [54]. Classical molecular dynamics simulations are based on the statistical mechanics. It is meaningless to focus on certain trajectory, and people are usually interested in the average values of a thermal dynamic system, such as energy, temperature and pressure. According to the different simulation set up, three ensembles are commonly used in the classical molecular dynamics simulations: (1) microcanonical ensemble (NVE), (2) canonical ensemble (NVT), and (3) isothermal-isobaric ensemble (NPT). Microcanonical ensemble denotes a system isolated with external environment, which has conserved atom/molecular number (N), volume (V) and energy (E). In canonical ensemble, the number (N), volume (V) and temperature (T) are conserved. The system energy is exchanged within a thermostat bath [43]. The popular techniques to control temperature in NVT ensemble include velocity rescaling, Berendsen thermostat [55], Nose-Hoover thermostat [56, 57], and Langevin dynamics [58]. Isothermal-isobaric ensemble represents a system with conserved number (N), volume (V) and pressure (P). 5 The pressure is controlled by a barostat bath, a similar concept with the thermostat bath. The widely used barostat baths are Berendsen barostat [55], Andersen barostat [59] and Parrinello-Rahman barostat [60]. 1.3. Outline of the thesis In this thesis, we will investigate the fracture and wrinkling in graphene with topological defects using atomistic simulations and continuum models. This thesis consists of six chapters and is organized as follows. Chapter 1 is an overview introduction about the mechanics properties of graphene with grain boundaries and dislocations. The methodology of molecular dynamics simulation and our research goal are also presented. In Chapter 2, we review the flaw tolerance/insensitivity phenomenon in nanoscale fracture, which has been observed in a wide range of materials. We will then investigate whether nanocrystalline graphene will have flaw insensitive fracture using large scale MD simulations. We will also study the effective fracture energy for the nanocrystalline graphene by combining MD simulations and continuum models. Chapter 3 focuses on the interaction between crack and typical topological defects, such as dislocation and grain boundary. We will set up a large graphene strip containing an edge crack and dislocation to investigate the interaction between them. The results based on MD simulations will be compared with the existing theoretical models, such as dislocation interacting with linear elastic crack and the Dugdale cohesive crack model. We will also carry out MD simulations for crack interacting with grain boundaries to estimate the toughness of grain boundary and predict the crack propagation path. 6 In Chapter 4, we will employ large-scale atomistic simulations and continuum modeling to analyze the defects controlled wrinkling in graphene. A benchmark for the deformation configuration and stress state of a single dislocation in a large grpahene will be first studied to show the pronounced 3D deformation and verify the continuum model proposed in the thesis. We will further design the graphene configurations under specific defect distributions such as those leading to a sinusoidal surface ruga1 or a catenoid funnel with both atomistic simulations and continuum models. In Chapter 5, we first give a brief introduction of the studies of curved graphene and phase field crystal (PFC) method. It will be shown that the PFC method can be used to generate initial atomistic structure of graphene conforming on curved surface. We will review the numerical schemes for solving the governing equation of PFC method and implement the finite element method (FEM) in the open source software, which is more flexible to deal with complicated boundary conditions. Preliminary studies for the fracture of a sinusoidal graphene sheet generating with PFC are also presented. Chapter 6 concludes the entire thesis with major scientific findings and future directions of the study. 1 The Latin word ruga is used to refer a large-amplitude state of a wrinkle, crease, ridge or fold [61]. 7 Chapter 2 Flaw Insensitive Fracture in Nanocrystalline Graphene 2.1. Introduction Chemical vapor deposition (CVD) method [5, 10] is still the current popular method to fabricate the large scale graphene, which will unavoidably create topological defects, such as dislocations and grain boundaries [14]. Damage of the graphene sample can also be induced by the process of transform and devices fabrication [17, 21]. On the other hand, extrinsic defects, like holes, have been introduced to tailor the functionality of the graphene, such as the electronic properties [62] and band structure [63], water desalination [64] and DNA sequencing [65, 66]. Therefore, it is crucial to understand the collaborative effects of different defects on the mechanical properties of graphene, like modulus and strength. Fracture is among the most prominent concerns about the mechanical behaviors of single and nano-crystalline graphene. It is known that the fracture strength of a material generally increases when the characteristic size (e.g., sample size, grain size) of the material is reduced [67]. One interesting phenomenon is that material failure may become insensitive to the presence of cracks when the sample size is smaller than a critical value [31]. A previous study has shown that the fracture of perfect graphene remains sensitive to the cracks or holes even down to the atomic scale [28]. An interesting open question is whether nanocrystalline graphene (nc-graphene), which contains a network of nanoscale 8 internal grain boundaries, could be made flaw tolerant and if so, what would be the critical condition for this to happen? This question is non-trivial due to the out-of-plane membrane deflections of graphene during deformation and fracture. In this chapter, we will first review the concepts of flaw tolerance/insensitivity and related studies in this area. We then establish a molecular model for the uniaxial tensile simulations of a nanocrystalline graphene strip containing a pre-existing hole. The elastic modulus, failure strength and flaw insensitivity of a nanocrystalline graphene nanostrip will be investigated, followed by some discussion and conclusion remarks at the end of the chapter. 2.2. Concepts of flaw tolerance/insensitivity Inspired by natural materials with superior strength such as bone, tooth and nacre, which share a common staggered nanoscale mineral-protein composite structure at the bottom level, Gao and collaborators [31] showed that the nanocomposites in nature exhibit a generic mechanical structure in which the mineral particles of nanometer size are selected to ensure optimum strength and maximum tolerance of flaws (robustness). To illustrate this concept, consider a cracked mineral plate with a surface crack (Fig. 2.1a), for which the failure strength  mf can be calculated from Griffith’s criterion of fracture as   mf  Em  ,  , (2.1) Em h 9 where E m is the Young’s modulus of mineral,  is the fracture surface energy and h is the thickness of the mineral plate. The failure strength for a defect free crystal is defined as the theoretical strength S . As shown in Fig. 2.1b, there exists a critical length scale Em h*   2 , (2.2) S2 below which the fracture strength of a cracked crystal is equal to that of a perfect crystal. This phenomenon, referred to as flaw tolerance, has now been reported in various materials, for example ceramics [68], biological composites [31], nanocrystalline metallic thin film [69, 70], amorphous carbon [71], filament networks [72], spider silk [73, 74], metallic glass[75, 76] and Pt nanocrystalline nanopillar [77]. Figure 2.1. Theoretical modeling of flaw tolerance. (a) Schematic drawing of a surface crack in a mineral plate. (b) Comparison of the fracture strength of a cracked mineral platelet predicted from the Griffith criterion of fracture with the strength of a perfect crystal. (a)-(b) are both from reference [31]. In general, the critical size for flaw tolerance depends on the sample geometry and loading conditions. Gao and Chen [78] conducted a more rigorous theoretical derivation of the critical size for flaw tolerance in a thin strip under tension, which provides a 10 benchmark solution for the verification of experiments or numerical simulations. The theoretical modeling set up is shown in Fig. 2.2, including a thin strip with a central crack (Fig. 2.2a) or an edge crack (F.2.2b). They defined the state of flaw tolerance as such that a pre-existing crack does not propagate even as the applied stress approaches the limiting strength of the material, denoted as S. There is no stress concentration in the structure so that the material failure occurs by uniform rupture at the limiting strength of material. Since the introduction of a crack will reduce the effective cross-section of material near the crack, the value of the far field applied stress should decrease as  a  th  S 1    S 1    , (2.3)  H where  th represents the theoretical failure strength in the presence of the flaw. Combining Eqs. (2.1) and (2.3) with considerations of the strip geometry, Gao and Chen obtain the critical size for flaw tolerance as H *  3.58 E S 2 , where   2 is the fracture energy of the material. 11 a c b Figure 2.2. Theoretical modeling of flaw tolerance in a cracked strip. (a) A thin strip with a central crack. The strip width is 2H and the crack length is 2a. (b) A thin strip with edge cracks. The strip width is 2H and the crack length is a. (c) Theoretical analysis of flaw tolerance in the thin strip. (a)-(c) are all from reference [78]. 2.3. MD Simulations of fracture in nanocrystalline graphene Our MD simulations of nc-graphene are performed using LAMMPS [53]. The interatomic force is described by the adaptive intermolecular reactive empirical bond order (AIREBO) potential [46]. It should be noted that the smaller cutoff distance in the switching function of AIREBO potential should be taken in the range of 1.92~2.0 Angstrom to avoid a known non-physical post-hardening behavior [15, 16, 54] while remaining consistent with DFT calculations [20]. For the present study, this cutoff 12 distance was set at 2.0 Å following a previous study [54]. The simulation sample is a nc- graphene nanostrip with the length of 60 nm, width of 20 nm and average grain size of 2 nm. To test the flaw sensitivity of the nc-graphene nanostrip, a circular hole with radius ranging from 1 nm to 6 nm is created at the center of the strip. Periodic boundary condition is imposed in the longitudinal (x) direction of the strip, while the other two directions (y and z) are kept free. During simulations, NVT ensemble is adopted to maintain the temperature at 300 K using the Nose-Hoover thermostat [56]. The tensile fracture behavior of the nc-graphene nanostrip is investigated under a constant stretching strain rate of 5×108 /s, while the Virial theorem is employed to calculate stress distribution in the sample. Figure 1a shows a typical atomic configuration of the nc- graphene nanostrip after equilibration at room temperature, where numerous wrinkles enhanced by atomic mismatches across grain boundaries can be clearly observed. These out-of-plane deflections are known to help minimizing the energy of the system and thermally stabilizing the membrane-like crystalline structures [32, 33]. Figure 2.3b-d illustrates various types of defects along the grain boundaries, such as the pentagon- heptagon pairs, polygons and vacancies. 13 Figure 2.3. Simulation sample of a nanocrystalline graphene nanostrip. (a) Contours of out-of-plane displacement after equilibration at 300 K. (b)-(d) Representative defect structures revealed by differences in potential energy, which include (b) pentagon- heptagon and pentagon-octagon pairs, (c) vacancy and (d) polygon. (Scale bar: 5 nm) Figure 2.4a depicts the tensile stress-strain curve of the nc-graphene nanostrip containing a circular hole of 1 nm radius. Three different regimes of deformation behaviors can be identified. In the first regime, the slope of the stress-strain curve increases as the applied strain increases, which is characteristic of entropic elastic behaviors in a thin membrane. Such entropic elasticity is attributed to wrinkling of nc- graphene due to the presence of defects and thermal fluctuation. Initially, there exist large wrinkles with RMS out-of-plane displacement of 0.83 nm in nc-graphene, which is substantially larger than the magnitude of finite-sized intrinsic ripples in the suspended single-crystalline graphene sheet [79, 80]. As the applied strain is increased, the nc- graphene gradually flattens under stretching, and the resistance to extension rises due to 14 the reduced degrees of freedom. In the second regime, the wrinkles are largely flattened by the applied stretching and the sp2 C-C bonds are being directly stretched, as evidenced by the fact that the stress becomes linearly dependent on the strain. Fitting the stress- strain curve gives an effective Young’s modulus around 432.26±10.57 GPa. In the third regime, the nc-graphene nanostrip fails due to fracture. Figure 2.4b-d shows a sequence of snapshots capturing crack initiation and propagation in nanocrystalline graphene. Typically, a crack initiates at a grain boundary and then grows along the grain boundaries or deflects into a grain. The sample eventually ruptures as the crack runs through the whole strip with little or no plastic deformation and the crack surface is nearly perpendicular to the loading direction, indicating brittle fracture. Notably, the fracture did not initiate at the hole but rather at some distance away from it, contradicting to the expectations from the classical theory of stress concentration at the hole boundary. We also performed MD simulations of tensile fracture in a nc-graphene nanostrip containing an edge notch (See Supplementary Materials), with similar observations that fracture occurred away from the notch. While these results are essentially reminiscent of recent experimental observations of flaw-insensitive fracture in nanocrystalline aluminum strips [69], multiscale simulations in model materials [81] and atomistic simulations in silk protein nanocrystals [74], this is the first time that flaw insensitive fracture has ever been demonstrated in high strength (tens of GPa) materials via atomistic simulations. 15 a b c d e Figure 2.4. Deformation and fracture behaviors of a nanocrystalline graphene nanostrip I with a center hole under longitudinal tension. (a) Typical tensile stress-strain curve. (b-d) 16 I I Tensile stress contours (unit: GPa) at ε=8%, 9.1%, 9.25% and 9.35%, respectively. In (b) and (c), micro-cracks (marked by squares) initiate from topological defects on grain boundaries. In (d) and (e), micro-cracks coalesce to form a big crack, eventually leading to failure. Note that the fracture surface is nearly perpendicular to the loading direction. (Scale bar: 5 nm). To understand the physical mechanism of flaw insensitive fracture in nanocrystalline graphene, we have also investigated the tensile fracture behavior of a single-crystalline graphene nanostrip with the same dimensions as our nc-graphene sample but without grain boundaries. Based on the crystal orientation of the graphene edges normal to the tensile direction, the simulated sample is classified into armchair and zigzag types. Figure 2.5a shows the relationship between the strength (taking as the peak tensile stress before fracture) of the single-crystalline graphene strip and the radius of the center hole. The results are consistent with those from multiscale calculations coupling quantum, molecular and continuum methods [28], as well as MD simulations on single-crystalline graphene with a center slit [82]. According to the classical Griffith model, the facture strength of a center-cracked strip is [83] 1 E  , (2.4) F   a where 2a is the crack size, E is the Young’s modulus,   2a / W ,  is the fracture surface energy and F   is a function of  reflecting the boundary effect,  F    1  0.025 2  0.06 4  sec 2 . (2.5) Taking the Young’s modulus and fracture surface energy of graphene to be 1 TPa and 10 J/m2, respectively, the predictions from the Griffith model are in excellent agreement 17 with our MD simulation results, as shown in Figure 2.5a. The value of fracture energy used here is close to that obtained from coupled quantum-atomistic simulations, as well as from classic MD simulations [28], although a drastically different value (0.1 J/m2) has also been reported from a molecular statics simulation [84]. Next we make a quantitative comparison between the tensile fracture behavior of a nc-graphene strip and that of a single-crystalline graphene strip. For every selected value of the radius of the hole in the nc-graphene strip, we perform independent simulations on ten randomly generated samples with different initial grain configurations but the same mean grain size from the Voronoi construction. The fracture strength is then taken as the average of all ten simulations to eliminate statistical fluctuations from different simulations. Figure 2.5b shows the variation of fracture strength with the hole radius. As the hole radius decreases, the fracture strength increases and then saturates to a plateau which is close to the theoretical strength (about 26.86 GPa) of nc-graphene without a hole. It is noted that a majority of previous studies focused on bi-crystalline graphene with specific arrangement of pentagons and heptagons along grain boundary [15, 16]. Due to the presence of relatively ordered pentagonal-heptagonal rings, the strength of such bi-crystalline graphene is higher than that of polycrystalline graphene with multiple grain boundaries of random misorientations. It is found that Eq. (2.4) is still capable of predicting the fracture strength of the nc-graphene strip if the fracture energy is reduced to 8 J/m2, and there exists a transition from the Griffith-governed fracture to failure at the limiting strength of material, as shown in Figure 2.5b. Such transition has also been demonstrated by numerical simulations for a thin strip containing central crack with lattice model [81] and a silk protein nanocrystals under point load [74]. Note that this is 18 the first time that the facture surface energy of nc-graphene is estimated. Due to the presence of grain boundaries, it can be expected that the fracture energy of nc-graphene should be a little smaller than that of pristine graphene. Figure 2.5. Strength variations of single- and nano-crystalline graphene nanostrips with the radius of an internal hole. (a) Strength of armchair and zigzag single-crystalline graphene.  EX is the experimental value of graphene strength reported in Ref. [4]. The reader is referred to Ref. [28] for strength calculated from a multiscale method coupling quantum, molecular and continuum statics and to Ref. [82] for strength from MD simulations of a square graphene sheet (55 nm2) with a central slit at room temperature 19 and strain rate of 5×108 /s. (b) Strength of nanocrystalline graphene. (c) Normalized strength of single- and nano-crystalline graphene nanostrips.  m denotes the strength from MD simulations,  th refers to the limiting strength of the flaw tolerant nanocrystalline graphene strip at the holed cross section. (d) Strengths of single- and nano-crystalline graphene strips with an ellipitic hole. According to the previous theoretical model [78], the critical width for a strip with a central crack to achieve flaw tolerance is Wcr  3.58 E S 2 , (2.6) where S is the limiting strength of material (i.e. the average strength of nc-graphene without holes). From the present MD simulations, the limiting strength of nc-graphene is obtained as 26.86±1.12 GPa. Therefore, the critical strip width for flaw-tolerance is estimated as 17.16 nm from Eq. (2.6). The width of the nc-graphene strips used in the present simulations is 20 nm, which is close to the critical value. When the strip width falls below the critical value, the strip becomes flaw tolerant and should fail at the limiting strength of material. To check this prediction, the strength is normalized as follows    m a   th a  , (2.7) where  m a  is the strength from MD simulation and  th a   S 1    is the limiting strength of the strip at the holed cross section. The relevant results for single- and nano- crystalline graphene strips are shown in Figure 2.5c. It is clearly seen that the normalized fracture strength of nc-graphene is independent of the flaw size, while the normalized strength of single-crystalline graphene falls way below the limiting strength. We have 20 also investigated the mechanical responses of single- and nano-crystalline graphene containing an elliptic hole of various aspect ratios, with results showing that the strength of single-crystalline graphene is sensitive to the aspect ratio of the elliptical hole, as would be expected from the classical theory of stress concentration. In contrast, the strength of the nc-graphene strip is found to be independent of the aspect ratio of the elliptical hole, again indicating flaw insensitive fracture. These results are summarized in Figures 2.5c,d, which highlight the differences in fracture behaviors between nano- and single-crystalline graphene nanostrips. Figure 2.6. Tensile stress contours (unit: GPa) near a hole of radius 3 nm in single- and nano-crystalline graphene nanostrips during deformation. (a) Armchair and (b) zigzag 21 graphene strips at a strain of 4%. Nano-crystalline graphene strips (c) with and (d) without the hole at a strain of 8%. (Scale bar: 5nm) Figure 2.6 compares the stress contours between single- and nano-crystalline graphene during deformation, showing clear stress concentrations near the hole in single- crystalline graphene. In contrast, while some elevated stress regions can be observed along the grain boundaries in nc-graphene, there is no evidence of stress concentration near the hole. This is consistent with recent experimental results on notched nanocrystalline aluminum thin strips [69, 70]. Compared with single crystalline graphene, the grain boundaries in nanocrystalline graphene increase the characteristic length scale of flaw tolerance while weakening the elastic modulus, fracture energy and strength of the material. Our findings are essentially consistent with previous studies demonstrating that hierarchical materials with weak internal boundaries can enhance the flaw tolerance of materials [31, 74]. Based on the present study and recent experimental observations [18] that cracks may penetrate grain boundaries and move along armchair or zigzag paths in polycrystalline graphene with mean grain sizes more than one hundred nanometers, it is speculated that polycrystalline graphene with larger grain sizes may have smaller flaw tolerance length scale compared to nanocrystalline graphene. The phenomenon of flaw insensitive fracture is not limited to the geometry of the flaw. To demonstrate this point, the same sample configuration depicted in Fig 2.4 is employed to investigate the tensile fracture behavior of a pre-edge-notched nc-graphene strip (see Fig. 2.7). The simulation settings are kept the same as in the case of a pre- existing hole. Figure 2.7 demonstrates that fracture could occur at a distance away from the pre-existing notch, similar to the case of a central hole discussed in section 2.3. 22 Figure 2.7. Flaw insensitive fracture in an edge-notched nc-graphene strip. (a) Typical tensile stress-strain curve of a pre-edge-notched nc-graphene strip. (b) Tensile stress contours (unit: GPa) of the strip at different strains. 23 2.4. Summary In summary, we have performed MD simulations to investigate the uniaxial tensile fracture behaviors of a nanocrystalline graphene nanostrip, with focus on the phenomenon of flaw insensitive fracture below a critical strip width. Our simulations reveal that micro-cracks nucleate randomly at intrinsic defects along grain boundaries and coalesce to form a big crack, eventually leading to catastrophic fracture. The key finding is that the microcrack nucleation and coalescence are not always induced by nor associated with the large hole or notch in the system. It is found that the MD simulation results are in good agreement with theoretical predictions based on fracture mechanics and also qualitatively consistent with experimental observations on nanocrystalline aluminum thin strips. The present studies also showed the entropic elastic behavior in the initial stage of deformation, and predicted the facture surface energy of about 8 J/m2 for nanocrystalline graphene, which is lower than that (about 10 J/m2) of single-crystalline graphene. These differences can be attributed to the presence of defects along grain boundaries in nanocrystalline graphene. It is interesting and worth noting that, while the grain boundaries substantially weaken the strength of nc-graphene, they simultaneously render the material less sensitive to structural flaws. Our simulations provide significant insights into previous experiments on flaw insensitive fracture in nanoscale materials and reveal a fundamental difference in tensile fracture behaviors between nano- and single- crystalline graphene. 24 Chapter 3 Crack Interacting with Dislocations and Grain Boundaries in graphene 3.1. Introduction In the previous chapter, the fracture behaviors of a nanocrystalline graphene nanostrip with and without holes have been studied. It is shown that grain boundaries can weaken the strength of the graphene, while enhance the flaw insensitivity due to the presence of intrinsic defects. However, the complicated grain geometry and random distribution of grain boundary orientation as well as many other types of defects (i.e. vacancies) make it hard to obtain a fundamental understanding of the toughening effects of dislocations and grain boundaries. Dislocation core energy is a crucial parameter to establish a consistent continuum model for dislocations [85] and can only be obtained from atomic simulations. So far, there are no systematic studies of dislocation core energy in graphene. It is believed that the pre-stress induced by dislocations is the key to determining the strength of a grain boundary [15] in graphene. The grain boundary toughness, another important fracture property of material, has not been discussed in the literature. A recent experiment measured the toughness of polycrystalline graphene via direct tensile loading and reported a toughness value of about 15.9 J/m2 [86]. A natural question is whether grain 25 boundaries weaken or enhance the toughness of graphene? Here we analyze some simple scenarios involving a crack interacting with a dislocation or a grain boundary, and hope these studies may shed some light on the inter- and intra-granular crack paths found in polycrystalline graphene [18]. This chapter is aimed to investigate the following questions through atomistic simulations: (1) what is the dislocation core energy for pentagon-heptagon pair in graphene? (2) Will a dislocation have a significant shielding effect on crack propagation in graphene? (3) What is the toughness of a grain boundary with different misorientation angles in graphene? 3.2. Dislocation core in graphene A dislocation dipole model is adopted here to compute the dislocation core energy, as shown in Fig. 3.1 a. Periodical boundary conditions are applied to a square graphene sheet. The dipole distance varies from b to 31b, where b is the Burger’s vector for a pentagon-heptagon pair in graphene. The conjugate gradient (CG) algorithm is employed to find the energy minimum configuration. The dipole energy is defined as the energy change due to the introduction of the dipole compared to the perfect graphene. In order to compare with the classical dislocation theory, the carbon atoms are only allowed to move in the same plane during the simulations (i.e., two-dimensional (2D) simulations). According to the classical dislocation theory, the strain energy stored in an infinite domain with a dislocation dipole under generalized plane stress condition can be expressed as [87], Sb 2 d Sb 2 E d   2 Ecore  n log  2 Ecore  log , (3.1) 4 rc 4  26 with the dipole distance |d|=nb and dislocation radius rc=b. In reality, the simulation is conducted within a large graphene supercell (20×20 nm2) to represent an infinite domain. L a b L d c d -7.8 -6.8 eV Figure 3.1 Dislocations dipole in graphene. (a), Molecular mechanics (MM) model for a dislocation dipole in a graphene square sheet (length L of 20 or 40 nm) with dipole distance d varying from 0.24 to 6 nm. (b), Potential energy for the dislocation dipole from MM simulations and classical edge dislocation model. (c), Atomic configuration of the dislocation dipole in graphene with d=26b. (d), Core of a heptagon-pentagon dislocation in graphene shown in (c). (Scale bar: 2 nm) 27 Fitting the MM simulation results for a dislocation dipole in the graphene supercell (see Fig. 3.1 b), we can obtain the following key parameters for the model, i.e., the dislocation core energy Ecore=2.386 eV for rc=b and a core radius of rc=0.59b=1.43 Å if we set Ecore=0 eV, which are in good agreement with previous values reported in literatures [88]. The total potential energy stored in the dislocation core (the region circled in Fig. 3.1d) is 2.362 eV, which is a self-consistent proof of the dislocation core energy and radius estimated from Eq. 3.1. The in-plane modulus S predicted by Eq. 3.1 is equal to 19.14 eV/Å2, very close to the value 17.57 eV/Å2 obtained from the uniaxial tensile simulation. The convergence of results with respect to the sample size is further confirmed by comparing the dislocation dipole energy in a larger supercell size (40×40 nm2 shown as red square in Fig. 3.1 b). The agreement between simulations and theoretical predictions shows that the classical dislocation model can be applied to describe dislocations in graphene confined to 2D. 3.3. Crack interaction with a dislocation in graphene Crack-like defects commonly exist during large area graphene fabrication and transformation processes. A fundamental question is how the interaction between cracks and dislocations affect the strength of polycrystalline graphene. Here we use MD simulations to investigate the interaction between a single dislocation and a crack. It is known that the relative orientation of a dislocation with respect to a crack tip plays an important role. Two typical dislocation orientations are considered here, i.e. dislocation-1 denotes the case of its pentagon facing the crack tip and dislocation-2 represents the case of its heptagon facing the crack tip. 28 a b d e c 1 nm f g h i 1 nm -80 80 GPa Figure 3.2 Crack interaction with a dislocation in graphene. (a), Stress-strain curves from MD simulations for a cracked graphene strip with/without a dislocation ahead of the crack tip. (b), Stress contours (yy) for a graphene strip with crack only (εyy=0.027). (c), A zoom-in snapshot around the crack tip in (b). (d), Stress contours for a crack interacting with a dislocation-1 (εyy=0.035). (e), A zoom-in snapshot around the crack tip in (d). (f), Stress contours for a crack interacting with a dislocation-2 (εyy=0.019). (g), A 29 zoom-in snapshot around the crack tip in (f). (h), Stress contours for a crack interacting with a dislocation-2 (εyy=0.019). (i), A zoom-in snapshot around the crack tip in (h). The scale bars are 10 nm in (b), (d), (f) and (h). A graphene nanostrip with length of 80 nm and height of 40 nm, and containing a 20 nm long edge crack created on the middle plane is used in current simulations. The initial crack remains atomically sharp and along the zigzag plane of graphene. The dislocation is placed 1 nm away from the crack tip. A tensile loading with strain rate of 109/s is applied by prescribing a linear velocity profile at the beginning and then the vertical velocities of the upper and lower edges of the strip. During simulations, NVT ensemble is adopted to maintain the temperature at 10 K using the Nose-Hoover thermostat [56]. The failure strain for the cracked graphene without dislocation is 0.0277 (see Fig. 3.2a), from which the fracture surface energy can be extracted as Γ =13 J/m2 (i.e. Γ=Ehε2, E being the Young’s modulus, h the half width of the strip, ε the critical strain at crack propagation). For the case of crack interacting with dislocation-1, the final failure strain increases to 0.036 (see Fig. 3.2a), indicating a significant shielding effect, which can be also verified by comparing the stress field near the crack tip (shown in Fig. 3.2b-e). For the case of crack interacting with dislocation-2, a tensile pre-stress near the heptagon will enhance the tensile stress near the crack tip and assist the crack initiation, as shown in Fig. 3.2f-g. The critical strain for initial crack propagation is about 0.02. The advancing crack tip will break the weak bond associated with the heptagon-hexagon pair, resulting in an effectively blunted crack tip and hindering the catastrophic failure of graphene strip until a much larger critical strain of 0.034. 30 It is interesting to compare our current MD simulations with the existing theories about crack interacting with dislocations. According to linear elastic fracture mechanics, the shielding effect of a dislocation on a crack can be expressed in term of an effective stress intensity factor at the crack tip [89]  sin  Eb 2  K Ieff  K I0 1  , (3.2)  8 R    where KI0 is the stress intensity factor for perfect without dislocation, E is the Young’s modulus, Γ is the fracture energy, b is Burger’s vector and R is the distance between dislocation and the crack tip, ϕ is the angle of Burger’s vector (ϕ= π/2 for dislocation-1, and ϕ= -π/2 for dislocation-2). At the atomic scale, the crack-tip singularity can be eliminated by the nonlinear deformation of atomic bonds. Bhandakkar et al [90, 91] studied a dislocation interacting with a Dugdale cohesive crack and derived a closed form formula for the effective stress intensity factor as  L sin  Eb 2  K Ieff  K I0   , (3.3)  c 8 R    and K Ieff L L b c   L   1 2   sin   2 tan 1    ,  (3.4) c c    K I0   c    E where  c  is the length scale associated with the Dugdale cohesive zone [92] and 8  c2 L is the cohesive zone length in the presence of the dislocation. The unknowns K Ieff and L should be solved by combining Eq. (3.3) and Eq. (3.4). The detailed derivation can be found in papers [90, 91]. According to the model, the failure strain is linearly 31 propositional to the stress intensity factor, and therefore we can estimate the failure strain once we have the effective stress intensity factor. With the following parameters obtained from MD simulations of graphene,   2  13 J m 2 , R  1 nm, E  845 GPa, b  0.242 nm, c  120 GPa , (3.5) we can calculate the fracture strains for the cases of dislocation interacting with a linear elastic crack and a Dugdale cohesive crack, which are summarized in Table. 3.1. It is interesting to note that the linear elastic crack model overestimates the effect of dislocation, while the Dugdale cohesive crack underestimates the influence of dislocation. This does not seem unreasonable, as both of these models are approximations of the nonlinear cohesive behavior of atomic bonds near the crack tip. Of course, more sophisticated models of dislocation interacting with a cohesive crack can be developed to better calibrate the MD simulation results. It is also noted that the crack tip can break the bond shared by hexagon and heptagon when it is close enough to the crack tip. Under that situation, two crack tips will be formed, leading to an effective blunted crack and thus delaying the ultimate failure of the whole structure. In this sense, the failure behavior of graphene can be quite sensitive to the atomic structure near the crack tip. Table 3.1. Simulated failure strains of dislocation interacting with a linear elastic crack versus a Dugdale cohesive crack from MD simulations and continuum models. MD Results Linear Crack Cohesive Crack Dislocation-1 0.036 0.0385 0.0344 Dislocation-2 0.020 0.0169 0.022 Note: The simulated failure strain of a cracked graphene strip without a dislocation is 0.0277. 32 3.4. Crack along a grain boundary in graphene The theoretical strength of a material is an indicator of how strong of a defect free sample is. So far, graphene has been reported to have the highest strength (130 GPa [4]) among all materials. In comparison, the toughness of a material describes the critical condition for a crack to propagate in the material. Recent experiments found that the toughness of graphene is only 15.9 J/m2 [86], indicating that graphene is quite fragile. Here we will focus on the calculations of grain boundary toughness of graphene through MD simulations. Our simulations use a bi-crystal graphene strip with length of 80 nm and height of 40 nm, containing a symmetrical grain boundary in the middle (as shown in Fig. 3.3a). The initially atomically sharp crack lies along the grain boundary with length of about 20 nm. For the same grain boundary, there exists two different types of crack according to its propagation direction, namely a rightward crack (Fig. 3.3b) and a leftward crack (Fig. 3.3c). A mode I tensile loading with strain rate of 109/s is applied similarly to the previous simulations. During the simulations, NVT ensemble is again adopted to maintain the temperature at 10 K using the Nose-Hoover thermostat [56]. 33 a b c -7.8 -6.8 eV Figure 3.3 Cracking along a grain boundary in graphene. (a), Bi-crystal graphene strip with left crack. (b), A rightward crack along a zigzag grain boundary (θ=9.43). (c), A leftward crack along an armchair-tilted grain boundary (θ=9.43). The scale bars are 10 nm in (a) and 1 nm in (b) and (c). 34 a b -6.8 eV c -7.8 d e ε=0.020 ε=0.020 f g ε=0.029 ε=0.032 h i ε=0.030 ε=0.034 -80 80 GPa Figure 3.4 MD simulations of a crack propagating along an armchair-tilted grain boundary with misorientation angle θ =9.43. (a), Stress-strain curves for the rightward 35 and leftward cracks. (b) and (c), Crack propagation patterns. (d), (f) and (h), Snapshots for stress (σ22) distribution near the crack tip for the rightward crack under different applied strains. (e), (g) and (i), Snapshots for stress (σ22) distribution near the crack tip for the leftward crack under different applied strains. The scale bars are 10 nm in (b-c) and 1 nm in (d-i). The results for crack propagation along a typical armchair tilted grain boundary with misorientation angle θ =9.43 are shown in Fig. 3.4. As shown in Fig. 3.4a, the failure strain for the rightward and leftward cracks are 0.03 and 0.034, respectively, which are both higher than that of crack propagation along the armchair direction in a perfect graphene sheet (0.0277). This indicates that the grain boundary can enhance the toughness of graphene, even though the strength is lower. Our simulations also show that the crack will keep propagating along the grain boundary during the whole process (Fig. 3.4 b-c). Some details of crack initiation and propagation can be examined by comparing a series of snapshots of stress distributions near the crack tip (Fig. 3.4 d-i). For the rightward crack, atomic bonds near the crack will start to break as early as 0.02 applied strain. However, the crack tip will be trapped as it approaches a pentagon and the applied strain will continue to increase (Fig. 3.4 f) until the bond shared by the nearest hexagon and heptagon is broken (Fig. 3.4 h). Although the leftward crack is similarly hindered by the grain boundary, the details of crack propagation are are quite different. As shown in Fig. 3.4e, the maximum stress is not at the crack tip but near the heptagon for the leftward crack. The bond shared by the nearest hexagon and heptagon ruptures before the main crack starts to propagate (Fig. 3.4g), effectively creating a daughter crack which propagates towards to the mother crack to induce an overall crack propagation (Fig. 3.4i). 36 To understand the mechanism of crack propagation along the grain boundary, the stress distributions near the crack tip for the rightward and leftward cracks before any external loading is applied are compared in Fig. 3.5. It is clearly shown that the grain boundary causes an initial tensile stress at the rightward crack tip while an initial compressive stress at the leftward crack tip, which is contrary to a direct superposition of the stress field induced by dislocations in the grain boundary. It is known that the terminated dislocation wall can be treated as a disclination [93]. For the rightward crack, the dislocation wall forms a negative disclination, resulting in a tensile stress field (Fig. 3.4c). For the leftward crack, the dislocation wall forms a positive disclination, leading to a compressive stress field (Fig. 3.4d). While the simulated results can be qualitatively explained by this effective disclination model, the non-trivial geometry makes it difficult to obtain a quantitative analysis. a b c d x2 x2 s l l s … … x1 x1 -10 10 GPa 37 Figure 3.5 Initial stress (σ22) distributions near the crack tip for a crack along an armchair-tilted grain boundary. (a) MD simulations of a rightward crack. (b), MD simulation of a leftward crack. (c), Continuum models for (a). (d), Continuum model for (b). The scale bar is 1 nm in (a-b). a b c d e ε=0.28 ε=0.25 f g ε=0.29 ε=0.26 -80 80 GPa 38 Figure 3.6 MD simulations of a crack propagating along an armchair-tilted grain boundary with misorientation angle θ =21.8. (a), Stress-strain curves for the rightward and leftward crack propagation. (b) and (c), Crack propagation patterns for the rightward and leftward cracks. (d) and (f), Snapshots of stress (σ22) distributions near the crack tip for the rightward crack under different applied strains. (e) and (g), Snapshots of stress distributions near the crack tip for the leftward crack under different applied strains. The blue lines in (b) and (f) represent the lowest energy surface. The scale bars are 10 nm in (b-c) and 1 nm in (d-g). We also conducted similar simulations for an armchair-tilted grain boundary with misorientation angle of θ =21.8, as shown in Fig. 3.6. It can be seen that the failure strain for the rightward crack is comparable with that without the grain boundary, while the value for the leftward crack is smaller than that without the grain boundary. An interesting phenomenon is that the rightward crack tends to branch and deflect into the lower grain (Fig. 3.6b), while the leftward crack moves along the grain boundary during the entire process (Fig. 3.6c). Directional anisotropy in failure mode has been previously discussed for symmetrically tilt grain boundaries in copper nanocrystals [94]. The snapshots of crack tip field suggest that the branching of the rightward crack follows the lowest surface energy direction (Fig. 3f). In comparison, the leftward crack advances by breaking the bond shared by hexagon and heptagons (Fig. 3g). It is worth making a brief comparison between the current simulations and previous studies on crack propagating along grain boundary in metals. In metals, grain boundaries are usually regarded as weak planes and sources of dislocations emission. These two competing mechanisms can determine the intrinsic brittle or ductile behaviors of the grain 39 boundary, as described by Rice’s model [95]. The rapid development of computational resources have made it possible to model directional anisotropy associated with crack propagation along grain boundaries [94], lattice trapping effects [96] and impurity embrittlement of grain boundaries [97]. A recent study finds that disclination formed by two grain boundaries can heal an edge crack in nanocrystalline metals [98]. Our present simulations reveal that grain boundary can play a critical role in determining crack propagation behavior in a bi-crystal graphene. The pre-stress induced by dislocations in the grain boundary can have significant effect on the behavior of the crack tip. Apparently, there exist two competing mechanisms to determine whether the crack might branch away from an armchair-tilted grain boundary in graphene: crack propagation along a plane of lowest surface energy versus breaking the bond shared by hexagon and heptagon. It may be possible to develop a Rice-like continuum model to characterize the competing effect of these two factors in the future. 3.5. Summary In this chapter, we performed MM simulations to search for the minimum energy configuration of a dislocation dipole in a square graphene sheet under periodical boundary conditions. By calibrating the MM simulated data with a continuum model, we can extract the effective dislocation core energy and core radius. As a consistent check of our results, the total potential energy change of atoms inside the dislocation core is verified to be close to the fitted value for the dislocation core energy. The dislocation core radius obtained in our simulations is also in good agreement with the results from previous studies in the literature. 40 Our MD simulations on dislocation interaction with an edge crack in a 40 nm wide graphene strip reveal that a single dislocation can increase or decrease the effective stress intensity factor by 30% depending on the orientation of the dislocation with respect to the crack. It is also interesting to note that the MD simulation results lie between those predicted for a linear elastic crack and for a Dugdale cohesive crack. Therefore, it may be worth extending the current model for dislocation interacting with a Dugdale cohesive crack to a more realistic cohesive crack model. The simulations also show that a slight variation in the atomic structure near the crack tip can significantly change the failure strain. We also carried out MD simulations for crack propagation along an armchair symmetrically tilted grain boundary. Similar to a dislocation, grain boundaries can strongly effect the stress intensity factor at the crack tip. However, the pre-stress obtained from MD simulation has an opposite sign compared to that calculated from a direct superposition of the stress induced by dislocations in the grain boundary. This can be qualitatively explained from the effective dislcination of a terminated dislocation wall in a consistent continuum model, even though it has been difficult to obtain quantitative agreement. Two competing mechanisms are identified to determine the crack branching behavior along an armchair-tilted grain boundary in graphene: cracking along a plane with the lowest surface energy versus breaking of a bond shared by hexagons and heptagons in the grain boundary. In this regard, our current studies are still quite preliminary, and more work is needed in the future. For example, more systematic studies on crack propagation along grain boundaries with different misorientation angles may help achieve a complete understanding of grain boundary fracture in graphene. It will be 41 interesting to study crack propagation along grain boundaries under mixed mode loading conditions, which may be helpful for understanding the failure of polycrystalline graphene. It is noted that there is still a general lack of reliable interatomic potentials to capture the fracture behavior in graphene. While the fundamental mechanism of crack interaction with dislocations and grain boundaries via the pre-stress distributions at the crack tip may be insensitive to the detailed forms of the atomic potential, the exact values of the critical strains and the sequences of the bond breaking processes near the crack tip may be more dependent on the choice of the interatomic potential. Our studies are not only suitable for graphene, but also for other material systems like nanocrystalline diamond and cubic-BN crystal with dislocations or grain boundaries. 42 Chapter 4 Defects controlled wrinkling and topological design in graphene 4.1. Introduction As illustrated in the title, this thesis is dedicated to investigating two main problems in graphene with topological defects: fracture and wrinkling. In the previous two chapters, we have focused on the fracture in graphene with topological defects through MD simulations and theoretical analysis. Our studies have revealed that the defects such as dislocations and grain boundaries have significant influences on the strength and toughness of graphene. In the following two chapters, we will investigate the defects controlled wrinkling in graphene and the topological design of multifunctional graphene with defects. As a special class of solid membrane with two-dimensional (2D) crystalline structure, free standing graphene undergoes fully three-dimensional (3D) deformation to minimize its energy in the presence of topological defects as verified by recent high- resolution transmission electron microscopy (HRTEM) experiments [36, 37], density functional theory (DFT) calculations [32] and atomistic simulations [33-35]. As illustrated in Fig. 4.1a, a dislocation dipole in graphene can induce large wrinkles near the dislocation core, with out-of-plane displacement amplitude up to 3.3 Å. Figure 4.1b shows the potential energy of carbon atoms around the dislocation dipole, where atoms in the dislocation cores exhibit higher energy than those in the far field. In particular, three 43 atoms on the heptagon side have the highest energy. The 3D distributions of bond lengths around a dislocation and the corresponding 2D projection are presented in Fig. 4.1c and 4.1d, respectively. It can be clearly observed that the 2D projection significantly underestimates the length of covalent bonds around the dislocation core [36]. Interestingly, it has also been shown that the out-of-plane displacement in graphene can be used to tune the electrical [38-40] and mechanical properties of graphene, such as anisotropic friction [41] and tunable wettability through controllable crumpling [42] of graphene on substrate. These recent studies are calling for systematic investigations of topological defects induced wrinkling in graphene. Figure 4.1. Atomic configuration of a dislocation dipole in a graphene sheet of dimension 20×20 nm2 (only part of the region around the dislocation dipole is shown here). (a) A perspective view of deformation around the dislocation dipole. (b) Top view of the dislocation dipole. (c) Bond structures around the dislocation core in 3D. (d) Bond structures around the dislocation core in 2D projection. The color represents the scale of the out-of-plane displacement in (a) and potential energy in (b-d), respectively. (Scale bar: 1 nm) This chapter is organized as follows. Section 4.2 reviews the existing continuum model of defects in a flexible solid membrane. A mathematical analogy between 44 topological defects and incompatible growth metric will be established in Section 4.3. The numerical methods used in continuum and atomistic simulations are summarized in Section 4.4. Section 4.5 is dedicated to the fundamental problem of an isolated dislocation in graphene. Section 4.6 presents two examples of designing graphene-based nanostructures with distributed topological defects, specifically a sinusoidal graphene ruga and a catenoid graphene funnel. Finally, some concluding remarks are given in Section 4.7. 4.2. Continuum model for graphene with topological defects From a continuum perspective, topological defects in flexible membrane and shell structures have been investigated during the past several decades. Nelson and collaborators [99, 100] derived a generalized von Karman equation to describe coupling between topological defects, in-plane stress and out-of-plane deformation in a 2D crystalline membrane; they showed that the out-of-plane deformation can significantly reduce the magnitude of in-plane stresses generated by the defects. Zubov [101-103] investigated topological defects in plates and shells, and proposed a geometrical analogy to transform the problem of a thin shell with defects into its dual problem of a thin shell with external loading [103]. In this section, we briefly review the governing equations of topological defects in a flexible solid membrane by Seung and Nelson [100]. Our present study focuses on the topological defects in a mono-layer graphene sheet, reference configuration of which is taken to be a flat surface in the (x1, x2) plane. The deformation can be characterized by 45 the in-plane displacement components u1(X) and u2(X), as well as the out-of-plane displacement w(X). The in-plane strain tensor can be expressed as,  ij  1 ui, j  u j ,i  w,i w, j , i, j  1,2 , (4.1) 2 where ,i denotes partial derivative ,i xi . The mean and Gaussian curvatures can be written in terms of the out-of-plane deformation w as,   w   det w,ij  H   ,  G    , (4.2) 2 2  1  w 2  1   w   where  denotes the gradient operator and det  the determinant of a second rank tensor. For small w , Eq. (4.2) can be simplified as H   2 w,  G  det w,ij  , (4.3) where  2  is the Laplace operator. The total strain energy of the membrane is expressed as a summation of the stretching energy and bending energy, 1 Eh  v  1  E     ij ij   kk  kk dx1dx2     BH 2  BG G dx1dx2 , (4.4) 2 1 v  1  2v   2  where E is the Young’s modulus, v the Poisson’s ratio, B the bending stiffness, and BG the Gaussian stiffness. Here we note that the thickness of graphene is not well defined within the classical theory of plates and shells, as it may depend on the type of loading and geometry [104, 105]. In the present study, the effective thickness h of graphene (3.34 Å) is only used to define the in-plane stretching modulus, while the bending stiffness is obtained directly from atomistic simulations [106, 107] . 46 Disclinations characterize the rotational defects in a crystalline structure, whose direction is defined by a rotation axis also known as the Frank vector. For disclinations in graphene, the Frank vector is normal to the x1-x2 plane, so that the out-of-plane displacement remains single-valued. In the presence of N disclinations, the incompatibility condition should satisfy [100], eik e jl  u j ,i  ui , j    si r  ri  , 1  N (4.5) 2 ,lk i1 where eij is the 2D permutation tensor, and si(r-ri) represents the ith disclination at position ri with strength si. In the present study, dislocations are modeled as disclination dipoles [16, 93, 108], which makes the distributed dislocations not explicitly appear in the incompatibility equation, implying e jl ui , jl  0 [100]. Substituting this into Eq. (4.5) and considering the identity 1 ui, j  u j ,i    ij  1 w,i w, j , we obtain the incompatibility 2 2 condition in terms of in-plane strains and out-of-plane deformation as follows,   N eik e jl   ij  w,i w, j    si r  ri  . 1 (4.6)  2  ,lk i 1 Taking variation of the total energy in Eq. (4.4) with respect to the in-plane and out- of-plane displacements, we obtain the following equations, B 4 w  h ij w, j ,i , (4.7)  ij ,i  0 where 4 is a bi-harmonic operator and the stress tensor is calculated according to E Ev generalized Hooke’s law  ij   ij    , for i,j,k=1,2. 1 v 1  v 1  2v  kk ij 47 Introducing the Airy stress function , so that ij=eikejlh-1,lkor11= h-1,22, 22= h- ,11 and 12 -h-1,12, the equilibrium equations ij,i=0 are automatically satisfied, and 1 the strains are expressed as 1 v v  ij  eim e jn  ,nm   2  ij , (4.8) S S where S=Eh denotes the in-plane stretching modulus. Substituting Eq. (4.8) into the incompatibility equation in (4.6) and using Eq. (4.3), we obtain  N   4    S  G   si r  ri  . (4.9)  i 1  Note that h w  ij , j ,i  eik e jl ,lki w, j  eik e jl ,lk w, ji  eik e jl ,lk w, ji . (4.10) The governing equation for nonlinear coupling between in-plane stresses and out-of- plane deformation can be rewritten as, B 4 w  w, , (4.11) where [f, g]= eikejlf,lkg,ji =f,11g,22+ f,22g,11 2f,12g,12. Finally, we obtain the following generalized von Karman equation for defects in flexible membrane, B 4 w  w,    N . (4.12)  4    S  G   si r  ri   i 1  It is noted that Chen and Chrzan [109] calculated the self-energy of a periodical array of dislocation dipoles in a graphene sheet with and without out-of-plane deformation by modeling dislocations as topological constraints and performing energy minimization in the Fourier space. In spite of these impressive progresses on dislocations/disclinations in solid membranes, it remains a challenge to develop an efficient continuum model that can 48 accurately capture both the atomic-scale rippling near the defect core and large-scale wrinkles in the membrane. 4.3. Mathematical analogy between defects and inhomogeneous growth In linear elastic solids, dislocations can be treated as an eigenstrain field that can be analytically or numerically treated using Green’s function method [110]. In comparison, topological defects in graphene present a greater challenge due to the complex nonlinear coupling between in-plane and out-of-plane deformations. It has been difficult to obtain analytical solutions even for a single, isolated dislocation in a membrane. Similar problems with incompatible deformation significantly influencing the final morphology of the structures occur in the growth of biological tissue [111-114] and the swelling of a soft membranes [115, 116]. A number of theoretical frameworks have been proposed to treat incompatible deformation in these phenomena, such as non-Eulidean plate theory [117] and multiplicative decomposition of deformation gradient [118-120]. A comparison between the governing equations for defects and growth in thin plates would suggest that the defects can be represented as an equivalent eigenstrain field, an idea that has been widely used in micromechanics [110, 121]. The governing equation for the inhomogeneous growth in thin films [111, 112] can be written as, B 4 w  w,     4   S  G  g , (4.13) 49 where g  11, 22   22,11  212,12 is the incompatibility metric due to the in-plane growth g g g or swelling. It is noted that the two sets of equations (i.e., Eqs. (4.12) and (4.13)) are N identical if one sets g    s  r  r . Since the growth strain components, i 1 i i  11g ,  22g and  12 , are prescribed values and independent of each other, there are no unique g combinations of growth strain for a given incompatibility metric field. In the present study, we set 12g  0, 11g   22g   g , and obtain a Poisson’s equation for  g , N     si r  ri  , 2 g (4.14) i 1 which has the following solution in an infinite domain, N si  g   log r  ri  constant , (4.15) i 1 2 where the constant term is determined by boundary conditions. For general computational domain and boundary conditions, the eigenstrain field can be obtained by numerically solving the Poisson equation defined in Eq. (4.14). Therefore, we will use the following growth strain to simulate the topological defects in graphene, N si  12g  0,  11g   22g   log r  ri  constant . (4.16) i 1 2 To eliminate the singularity at the core of the defects, we replace  r  ri  with a 1  r  ri 2  Gaussian function exp   by introducing an intrinsic length scale rc, and rc2  rc 2  then modify the solution to Eq. (4.16) as, 50 si  1  r  ri   2   log  r  ri   Ei    constant , g (4.17) 2  2  rc2  where Ei(x) is the exponential integral. For discrete systems like graphene or finite element method discretized structures, the length scale rc can be physically related to the lattice length or the minimum mesh size. 4.4. Numerical methods for simulation To provide some benchmark solutions to validate the continuum model, we first investigate dislocations in graphene using atomistic simulations based on the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) [53]. The interatomic force for graphene is calculated by the Adaptive Intermolecular Reactive Empirical Bond Order (AIREBO) potential [46]. During simulations, energy minimization is employed to find the equilibrium configurations of the defective system. Before equilibration, the simulated system is relaxed at a constant temperature of about 1 K via the Nose−Hoover thermostat [56] to introduce some perturbations in the out-of-plane deformation. In the case of a single, isolated dislocation, a small perturbation is applied to the structure by initially assigning some non-zero out-of-plane coordinates for one or two atoms near the heptagon and then relaxing the whole system, as the thermal fluctuation breaks the initial flat configuration and structural symmetry. The atomic stress is calculated by the Virial theorem. The continuum model is numerically solved based on a discretization of graphene into a triangular lattice [100], in which the total elastic energy is defined as a combination of stretching energy Fs and bending energy Fb, 51 S  rij  r0  , Fb  B 1  n  n   , 3 2 F  Fs  Fb , Fs  2 (4.18) 4 ij 3  where rij is the current bond length, r0 the equilibrium bond length, and n and n the normal vectors of nearest-neighbors (Fig. 4.2a). It is worth noting that the above triangular model predicts that the Gaussian stiffness is equal to the negative bending stiffness [100]. This is close to a recent DFT calculation [122] which predicted the bending and Gaussian stiffnesses of graphene to be 1.44 eV and -1.52 eV, respectively. During simulations, the growth strain g is applied to the lattice by stretching the equilibrium bond length r0 to r0 [1+g(x1,x2)] [112]. The equilibrium bond length r0 is taken to be 1.40 Å, the same as that used in atomistic simulations based on the AIREBO potential. Figure 4.2b shows the predicted 3D buckled shapes of positive and negative disclinations in a circular disk from the continuum model, which are shown to be consistent with atomic simulations in Fig. 4.2c. Figure 4.2. Buckled shapes around disclinations in a circular graphene disk of radius 5 nm. (a) The triangular lattice model adopted in continuum simulations. (b) The calculated 3D buckled shapes around positive and negative disclinations in the circular graphene disk from the continuum model. (c) Corresponding 3D buckled shapes from atomistic 52 simulations. The colors represent eigenstrain field with rc=0.5 Å and out-of-plane deformation in (b) and (c), respectively. 4.5. An isolated dislocation in graphene In this section, both continuum model and atomistic simulations are used to study an isolated dislocation in a square graphene sheet with in-plane dimensions of 80×80 nm2, as schematically shown in Fig. 4.3a. In both atomistic and continuum simulations, the edge of graphene is clamped to reduce the edge effect on wrinkling. The stretching modulus and bending stiffness are extracted from atomistic simulations as S=17.57 eV/Å2 and B=0.95 eV, respectively, and then used in continuum simulations. According to the atomic structure of graphene, the strength and size of the disclination dipole that forms a dislocation are taken as π/3 and 1.158 Å, respectively. Figure 4.3b depicts the atomic structure of the heptagon-pentagon dipole, while Fig. 3c shows the eigenstrain field of the dipole used in the continuum model, which is calculated according to Eq. (16) with one positive and negative disclination at the position of (-1.158,0) and (1.158,0), respectively. Figure 4.3. An isolated dislocation in a square graphene sheet. (a) Schematic illustration. (b) A magnified view of atomic structure of a dislocation (heptagon-pentagon dipole) in graphene. (c) The eigenstrain field of dislocation used in the continuum model. (Scale bar: 1 nm) 53 Figures 4.4a-c and 4.4d-f plot the in-plane stress contours of an isolated dislocation from continuum modelling and atomistic simulations, respectively. Comparison between Figs. 4.4a-c and 4.4d-f indicate that the stress field predicted from the continuum model is in good agreement with that from atomistic simulations. In order to further quantify the influences of the out-of-plane deformation on the stress field around dislocation, we plot stresses predicted from 2D atomistic simulations with no out-of-plane deformation and from the classical plane-stress solution of an edge dislocation in Fig. 4.4g-i and 4.4j-l, respectively. The results from 2D atomistic simulations are very close to the plane-stress theoretical solutions of an edge dislocation [85], but are remarkably different from 3D simulation results shown in Figs. 4.4d-f. It is observed that the stress fields with and without the out-of-plane deformation deviate significantly in magnitude, shape and distribution. For example, when the out-of-plane relaxation is allowed, the magnitude of the maximum tensile (compressive) stress 22 is seen to drop from 117 (-109) GPa to 58.5 (-48.6) GPa. Compared to the 2D solution, the region of significantly amplified stress near dislocation core is significantly reduced in 3D due to out-of-plane rippling. These results indicate that the out-of-plane wrinkling in graphene substantially lowers the stress level near the dislocation core, and the stress field in 3D becomes asymmetric due to the buckling prone heptagon (compressive) side of the dislocation core. 54 Figure 4.4. Comparison of stress fields around an isolated dislocation calculated from continuum modelling and atomistic simulations. Contour plots of stress components 11,22, 12 from (a-c) continuum modelling based on the generalized von Karman equation; (d-f) 3D atomistic simulations; (g-i) 2D atomistic simulations (without out-of- 55 plane deformation); (j-l) classical plane-stress solution of an edge dislocation. (Scale bar: 1 nm.) Figure 4.5 shows excellent quantitative agreement in out-of-plane displacement profile calculated from the continuum model and atomistic simulations. It is interesting to observe that the out-of-plane rippling extends over large distances (more than 20 nm), almost to the clamped boundary. This long-range nature of the out-of-plane distortion induced by a dislocation is consistent with recent HRTEM observations [36, 37]. Figure 4.5. Comparison of the out-of-plane displacement field calculated from continuum modelling based on the generalized von Karman equation and atomistic simulations. (a) Contour plots of the out-of-plane displacement field calculated from continuum modelling. (b) A zoom-in snapshot of deformation around the dislocation core in (a). (c) Corresponding contour plots of the out-of-plane displacement field calculated from atomistic simulations. (d) A zoom-in snapshot of deformation around the dislocation core in (c). The scale bars are 10 nm in (a) and (c), 1 nm in (b) and (d), respectively. 56 Most existing theoretical models of dislocations in graphene assumed 2D deformation and neglected out-of-plane displacement [16, 108, 123-125]. It has been argued that a 2D disclination dipole model can be used to predict the tensile failure strength of polycrystalline graphene with specific grain boundaries, with results in good agreement with those from molecular dynamics (MD) simulations [16, 108]. The reason for such agreement has been attributed to the fact that the wrinkles induced by dislocations can be suppressed by the applied tensile strain [16]. Our present model allows one to quantify the effect of an applied tensile strain on the wrinkling field around a dislocation in graphene. Figure 4.6. Effect of an applied tensile strain on wrinkling around a dislocation. (a) Continuum solutions under uniaxial and biaxial strains. (b) Corresponding atomistic simulations. (Scale bar: 1 nm) 57 Three different loading modes (uniaxial tension in x1, uniaxial tension in x2 and equibiaxial tension) are applied to an 80×80 nm2 graphene sheet that contains a single dislocation. Figure 4.6 shows the results from continuum modelling and atomistic simulation under the three loading modes with a maximum strain of 0.02. It is noted that only biaxial tension can effectively suppress the out-of-plane wrinkling. It appears that uniaxial tension in the x1 direction even enhances the out-of-plane deformation due to Poisson contraction [126]. The local wrinkling near the heptagon core is still suppressed by the uniaxial tensile strain. These results indicate a need for caution when using 2D models to analyze the deformation of dislocations in graphene even under tensile loading. It can be noted that there exists long-range wrinkling perpendicular to the loading direction under uniaxial tension, as shown in Fig. 6. Such wrinkling does not significantly affect the local stresses and strains near the dislocation, which determines the fracture strength of defective graphene [16, 22]. In this sense, the wrinkling under uniaxial tension has very minor influence on the strength of graphene. We note that the continuum model developed in the present study only involves three parameters, stretching modulus, bending stiffness and bond length (or Burger’s vector), which can be obtained from experimental measurements [4] or computational studies, such as DFT [107, 122, 127] and molecular mechanics [105, 106], with accuracy comparable to that of an interatomic potential [46]. Compared to full atom simulations, the continuum model can dramatically reduce the computational cost. For instance, typically it takes about 1000 CPU hours to perform atomistic simulations for an 80×80 nm2 graphene sheet, but only 50 CPU hours to perform continuum simulations of 58 comparable scale. The continuum model provides an opportunity to investigate defects controlled wrinkling in graphene at the micron scale and above. 4.6. Defects controlled graphene structures The nonlinear coupling between topological defects and curvature in a 2D crystalline membrane can induce a variety of interesting phenomena. It has been shown that lattice defects tend to adopt specific patterns on curved surfaces [128-134] and can trigger different buckling modes in a spherical elastic shell, which has been used to shed light on the understanding of virus shape [135] and morphological changes between smooth and faceted structures [136, 137]. Reversible transformation between flat sheets and surfaces with non-zero Gaussian curvature in nematic glass sheet has been studied by combining the effects of disclinations with thermal/optical stimuli [138, 139]. In parallel with recent attempts to design 3D surface profiles of inhomogeneous gel composites with controllable swelling ratio [115, 116], it will also be challenging and exciting to explore the possibility of designing curved structures with distributed topological defects in graphene and other 2D crystalline materials. The rapid developments in new experimental techniques [25, 140] have made it increasingly possible to control the position, type and distribution of defects in graphene, so as to tune its electrical and mechanical properties. The question of how graphene deforms due to isolated and/or collective defects in 3D motivates the present study to develop an efficient continuum approach to characterize 3D deformation and stress fields in graphene with prescribed distributions of defects, and to demonstrate the possibility of designing graphene-based nanostructures by controllably deploying such defects. 59 In the above discussions, we have demonstrated that significant out-of-plane wrinkling deformation in graphene can be triggered by topological defects such as disclinations and dislocations. One question is whether the collective behaviors of such topological defects can be utilized to design curved graphene structures with interesting properties. As a preliminary study in this direction, here we show two examples of curved graphene structures that can be generated and tuned by a distribution of topological defects. The present study can be related to the problem of finding the minimum energy configuration of a lattice of charged particles on a curved surface [130, 134]. Using Monte Carlo simulations, Hexemer et al. [130] predicted the emergence of disclinations in a lattice of charge particles on a high aspect ratio (amplitude over wavelength) sinusoidal surface wx1 , x2   A sinqx1 sinqx2  , (4.19) where A is the amplitude and q is the wavenumber. Figure 4.7. A sinusoidal graphene ruga induced by a periodic array of disclinations from continuum modeling and atomistic simulations. (a) The continuum eigenstrain filed in the 60 undeformed configuration with rc=0.5 Å and (b) the configuration of the sinusoidal graphene ruga predicted from continuum modelling. (c) Potential energy contour (top view) and (d) the configuration of the sinusoidal graphene ruga from atomistic simulations. (Scale bar: 1 nm) A sinusoidal graphene, referred to as a graphene ruga, can be created from a periodical array of distributed disclination quadrupole in a square unit cell [130]. Here we first create the graphene ruga in the continuum framework based on such a distribution of disclinations, whose corresponding eigenstrain field is shown in Fig. 4.7a. The continuum model successfully produces a sinusoidal ruga structure with amplitude of about 7 Å and wavelength of around 30 Å, as shown in Fig. 4.7b. Next we create an atomic graphene ruga through voronoi diagram of the lattice structure presented in the work of Hexemer et al. [130] and relax the structure using atomistic simulations under periodical boundary conditions. The atomic potential energy contour of the graphene ruga structure is depicted in Fig. 4.7c, from which we can see that most of the atoms have lower energy than those in the heptagon rings. Figure 4.7d shows the atomic sinusoidal ruga structure, which is nearly the same as the continuum prediction and thus validates that the continuum model can capture large scale wrinkles induced by multiple defects. More recently, Kusumaatmaja and Wales [134] showed that an array of dislocations emerges when particles confined on a catenoid surface are interacting with different potentials (Fig. 4.8a). The catenoid surface is defined as, x, y, z   c cosh ccos , c cosh csin ,  , (4.20) where 0 ≤ ζ< 2π, -zmax ≤ ξ < zmax and c is a parameter to control the waist of the catenoid in the z=0 plane. Typical catenoid surfaces with different waist are shown in Fig. 61 8a. To form a catenoid graphene funnel, an initial atomic structure of 1018 carbon atoms is constructed through voronoi diagram according to the coordinates of 600 triangular vertices from the Cambridge Energy Landscape Database (http://www- wales.ch.cam.ac.uk/CCD.html), corresponding to a shape determined from Eq. (4.20) with c=0.6. Before performing energy minimization via atomistic simulations, the edge atoms are functionalized with hydrogen atoms for passivation. The relaxed atomic structure is shown in Fig. 4.8c, confirming that the nearly catenoid graphene funnel constructed from topological defects is stable. To mimic the dislocation structure in the atomic system, we distribute 7 dislocation pairs (each composed of 2 heptagon-pentagon dipoles) on the upper and lower half of a cylindrical surface, as illustrated in Fig. 4.8d. The final configuration from the continuum modeling is shown in Fig. 4.8e, which is consistent with the full-atom model constructed from the catenoid surface. Figure 4.8. A catenoid graphene funnel from atomic simulations and continuum modeling. (a) Catnoid surfaces for c=0.5 and 0.7, from which zmax is determined by 62 cosh(zmax/c)c=1. (b) The initial atomic configuration of 1018 carbon atoms on a catenoid surface (Kusumaatmaja and Wales [134]) with distributed dislocations. (c) The relaxed atomic structure of (b). The color represents the potential energy scale. The corresponding (d) initial and (e) final structures predicted from continuum modelling. The color represents the eigenstrain scale with rc=0.5 Å. (scale bar: 1 nm) 4.7. Summary We have developed an efficient continuum approach to determining the 3D configuration and stress field of a curved graphene membrane that contains a distribution of topological defects. Comparison with large-scale atomistic simulations shows that the proposed continuum model is capable of predicting the stress field and out-of-plane deformation around disclinations and dislocations in graphene. Both atomistic and continuum simulations indicate that, due to the atomic thickness of graphene, even a single disclination/dislocation can cause significant out-of-plane wrinkling. The driving force for such wrinkling comes from the relaxation of in-plane stress and elastic energy near the defect core. We further showed that the out-of-plane wrinkling near a dislocation core cannot be fully suppressed by uniaxial stretching, which calls for caution in applying 2D dislocation models to analyze deformation around topological defects in graphene even under severe stretching. It should be possible to extend the present continuum method to other 2D one-atom-thick materials, such as monolayer B-N sheet and MoS2 membranes. A promising future research direction is to investigate how to utilize the nonlinear coupling between topological defects and membrane curvature to tailor the design of a curved graphene membrane or structure that conforms to an arbitrary three dimensional 63 surface or object. This is a highly nonlinear inverse problem, for which it would be prohibitively expensive to use full-atom simulation methods. In contrast, our present study has demonstrated the feasibility of developing a continuum framework which can accurately capture both the atomic scale rippling near the defects core and the large-scale wrinkling induced by multiple defects. The use of dislocations/disclinations to control configuration and stress/strain fields provides an extra dimension in design space to tune the electrical, thermal and mechanical properties of 2D crystalline nanomaterials. A highly efficient and accurate continuum approach may open the possibility to address a series of fundamental problems with immediate applications in various fields. 64 Chapter 5 A Phase Field Crystal Method for Multifunctional Curved Graphene Design with Topological Defects 5.1. Introduction In chapter 4, we have seen that 3D curved graphene sheets can be formed with distributed topological defects. It is believed that both the curvature and topological defects can alter the electronic [11, 14, 40] and mechanical properties of graphene [15- 17]. Thus, it is of great scientific and technological importance to develop a general methodology for creating a curved graphene sheet conforming to an arbitrarily given 3D shape for multifunctional device design. Our proposed continuum model can accurately capture both the global wrinkle and atomic configuration for a given distribution of defects. However, the defect distribution that generates a specific 3D shape of graphene membrane is usually unknown, which is actually an inverse problem involving highly nonlinear deformation. Previous studies have employed geometry method [141-143] and Monte Carlo simulations [144] to search for the positions of carbon atoms on a curved surface, which are difficult to apply in a large system with complicated geometry. One possible way to deal with this issue is to incorporate dislocation nucleation and motion based on our proposed continuum model in the last chapter. 65 It is noted that the recently developed phase field crystal (PFC) method [145] is a powerful promising approach to searching for the minimum energy configuration of a triangular lattice network on an arbitrary 3D curved membrane. One distinct advantage of the PFC method is that particle interaction is described by a continuum effective density field governed by a generalized diffusion equation, the time scale of which is several orders higher than MD simulations based on a pair potential. For example, Voigt’s group has successfully applied PFC to solve the Thomson problem of charged particles on a sphere [146, 147] and defect pattern of particles on a catenoid surface [148]. The triangular lattice structure obtained from PFC can be easily transformed into the hexagon lattice structure of graphene through voronoi construction. We will therefore use the PFC method to generate graphene atomic structure for any given curved membrane shape. In this chapter, we will first review some previous studies on the mechanical and physical properties of curved graphene and the methods used to generate such graphene sheets. We will then introduce the phase field crystal model and its implementation within the finite element (FEM) formulations to numerically solve the partial differential equation (PDE) associated with the PFC method. Finally, we will present a few preliminary results of curved graphene design with the PFC method. 5.2. Previous studies on curved graphene Curved carbon nanostructures have been found much earlier than Graphene. A famous example is fullerene (C60) (Fig. 5.1a), containing 12 pentagons to form a closed spherical shape, the discovery of which initiated the development of nanotechnology. The defects configuration on these closed geometries satisfy the classical Euler’s formula, given by 66 V  E  F  2, (5.1) where F, E, V are the number of faces, edges and vertices, respectively. If we assume all the carbon-carbon are sp2 bonding, the relationship between the numbers of different types of topological defects can be expressed as [149] 2 N 4  N 5  N 7  2 N 8  12 , (5.2) where Ni is the number of carbon rings with i edges. People have also tried to explore more complicated carbon material with curvature beyond the sphere. For example, Terrones and Mackay [150] proposed a three- dimensional, open and negatively curved graphitic structure. Lenosky et al. [151] conducted molecular mechanics simulation based on force field to investigate the energy of similar carbon structures and the associated defect structures. The negative curvature graphitic structure was subsequently identified by Iijima et al. [152] with transmission electron microscopy. Petersen et al. [144] developed a Monte Carlo simulations technique to investigate the complicated carbon nanostructured surface based on the reactive atomic potential. The readers can find more information from a nice review about many previous studies on curved carbon materials [149]. Very recently, experimental and simulated studies [153, 154] have shown that such kind of structure can help design carbon based supercapacitors. 67 a b c Figure 5.1. Curved carbon nanostructures. (a) Atomic structure of fullerene (C60) (http://en.wikipedia.org/wiki/Fullerene). Two-fold symmetry axis views of the cubic D (b) and P (c) of schwarzite, emphasizing differences in the layout of the seven-membered rings [151]. Since 2004, Graphene, a single atomic layer membrane, has attracted tremendous attention due to its exceptional electrical and mechanical properties [1-4]. The two- dimensional membrane structure of graphene presents a great challenge for real applications of graphene based devices, as it is difficult to manipulate graphene. It would be desirable to create some three-dimensional nanoporous form of graphene, a bulk material that has similar electrical properties as graphene. It has been reported that such kind of nanoporous structure (shown in Fig. 5.2) can be successfully fabricated with chemical vapor deposition (CVD) on a Ni-based metal foam structure [155]. It is believed that the nanoporous graphene structure can make it simpler to design novel electronic devices. 68 a b c Figure 5.2. Three-dimensional nanoporous graphene. (a) Fabrication of nanoporous graphene with the use of nanoprous Ni. (b) The structure of nanoporous graphene with and without nanoporous Ni. Raman spectral indicates that the 3D nanoporous sample consists of high-quality monolayered graphene. (c) Atomic structures of the nanoporous graphene. Topological defects can be found in the curved surface [155]. a b c d e 69 Figure 5.3. Creating a curved carbon structure from a geometrical point of view. (a) Unit cell for a toroidal carbon nanotube (Ref. [141]). (b) Bead molecular structure for a toroidal tube. (c) Bead molecular structure for hierarchical fullerene. (b) and (c) from the online blog: http://thebeadedmolecules.blogspot.com/. (d) Constructing curved graphene sheet from the triangular lattice network via the voronoi scheme. (e) Complicated curved carbon structures based on the triangular lattice network. (d) and (e) are from Ref. [143] Inspired by the curved graphitic structure connected with sp2 carbon-carbon bond, Bih-Yaw Jin, a professor in the chemistry department of National Taiwan University proposed to design and create 3D curved graphitic model structures using beads. They tried to combine chemistry, geometry and art together with bead molecular models; the structures they created serve as vivid demonstrations of the atomic structures of carbon based nanomaterials to the students. Jin’s group [141, 142] also derived the generalized classification scheme of toroidal and helical carbon nanotubes (shown in Fig. 5.3a) as well as other curved carbon materials from the geometrical point of view. The bead molecular structures for toroidal carbon annotubes and hierarchical fullerene are shown in Fig. 5.3b-c, which can be found on http://thebeadedmolecules.blogspot.com/. Similar to Jin’s concept of using beads to mimic carbon atoms, others [143] also constructed carbon atomic structures employing the voronoi construction from a triangular mesh (as shown in Fig. 5.3.d-e). For a given curved structure, the authors of [143] first used the sphere packing method to efficiently obtain a high quality triangular mesh with constrains that all the voronoi polygon have edges between 4 and 8. The curved carbon structures obtained from this method are confirmed to be stable using full atom simulations. This provides a reliable and efficient way to generate curved carbon nanostructures. The key 70 for this method is to create the initial triangular lattice structures, a problem which we will use the phase field crystal method to solve in our studies. 5.3. Phase field crystal method Phase field crystal is a recently proposed method, still under development, to describe the dynamics of discrete atomic structure with a continuum density function with periodical solutions, which can be obtained through energy minimization of a free energy functional and solved by a diffusion equation [145, 156]. Compared to the conventional continuum models and MD simulations, the PFC method can not only capture the atomic level information in a material, such as point and line defects, but also describe the dynamical process in real time scale. The classical functional used for PFC has been generalized from the Swift-Hohenberg (SH) equation [145]: 1  2  1  F        1       4 dx , (5.3) 2 4    where    2 is the Laplace operator, ϕ is the density, ε is the reduced x 2 y temperature representing the effect of undercooling. The governing equation for the density evolution can be expressed as:  t         1       3 . 2 (5.4) As shown in Fig. (5.4a), the conventional phase field model can only capture the amplitude variation of two phases, while the PFC method can track atomic structures. It has been shown that the parameters ε and  (the average density) determine the final patterns of density distributions [145, 157, 158], such as those in liquid, strip and triangle pattern for two-dimensional systems, as 71 liquid   , strip  As sinqs x    , and  triangle  At cosqt x  cos qt y   3  cos 2qt y   3 2  , (5.5) where q and A are the wave number and amplitude corresponding to the minimum of the free energy, 4 1 . qs  1, As  2  3   , qt  3 2 , At    15  36  2 2  (5.6) 5 3  The free energy of the system associated with each single mode approximation can be expressed as [157, 158]   2 4 Fliquid  1     , 2 4 2 1     2  4 Fstrip    5 , (5.7) 6 2 4  13 4   4 2 4  2  Ftriangle  0.1  2     15  36    . 2  50  2 25  5 3   The corresponding phase diagram for the 2D system based on this single mode assumption has been analyzed [145, 157] and summarized in Fig. 5.4b. For the triangular pattern, the system can be equivalent to an isotropic elastic system with effective elastic modulus calculated as 2 1   , and v  1 3 ,  3   15  36  2 (5.8) 75   where  and v are shear modulus and Poisson ratio, respectively. As the PFC method can capture the motion of topological defects (i.e. dislocations and grain boundary) in large time scale (seconds), it has been widely used to simulate the crystal grain growth [145] and grain boundary motion [159-161]. 72 a b Figure 5.4. The phase field crystal method. (a) Differences and similarities between the phase field and phase field crystal methods (from Prof. Voorhees’s talk at NIST). (b) Phase diagram of the phase field crystal model (Ref. [158]). The original PFC method adopts a very simple model that can generate periodical patterns to capture elastic and plastic (associated with defects motion) deformation of a discrete atomic lattice point. Later, people realized that the PFC model is actually a simplification of the classical fluid density functional theory (F-DFT) by perturbing the functional formula in F-DEFT around the average density [156, 162]. This connection not only provides a solid foundation of the PFC model, but also a systematic way to generalize the original PFC model to incorporate more complicated situations, such as multicomponent interactions [163], accurate models for 3D FCC and BCC crystal model [164] and new periodical patterns, i.e. hexagon and kagome shape [165]. With these new developments, it is now believed that the PFC model will play increasingly important roles in the large time scale dislocation dynamics in a number of materials, including FCC and BCC metals, and graphene [166]. In a real application, it is a still challenge to numerically solve the governing equation of the PFC model defined in Eq. (5.4), which is a 6th order nonlinear partial 73 differential equation (PDE). The higher spatial derivative not only requires special treatment in discretization methods in space, but also needs implicit solver of the time integral, as a sufficiently small time step must be used to make the integral stable (roughly t  x ). Currently, spectral methods in Fourier space are still popular in 6 solving PFC problems under periodical boundary conditions. However, it is known that numerical methods in real space, such as the finite element method (FEM), finite difference method and finite volume method, can be more flexible in handling the complicated geometrical shapes in solving general PDEs. Numerical schemes based on a mixed form of FEM, introducing two new variables to reduce the highest order in spatial derivative, have been proposed to solve the PFC equations [167], such as    , t   2u  u     3 , (5.9) u   , where μ and u are two newly introduced intermediate variables. If we are interested in the dynamical process on a curved surface  in 3D, like the defects in graphene considered in the current thesis, Eq. (5.9) can be rewritten as [147]    , t   2u    u     3 , (5.10) u    , where   is the surface Laplace operator. It has been shown that FEM can be used to successfully solve Eq. (5.9) and Eq. (5.10) [147, 158, 167]. Typical simulations for the PFC model on curved surface are shown in Fig 5.5. With FEM solver, Backofen et al.[147] applied the PFC method to the classical Thomson 74 problem of charged particles interacting on a sphere (Fig. 5.5 a) and obtained results consistent with previous studies in the literature. In the same paper, they also computed the lattice pattern on the Stanford bunny (Fig. 5.5 b) to demonstrate the advantage of the FEM scheme. Recently, Schmid and Voigt [148] analyzed the defect patterns on the catenoid surface with the PFC model (Fig. 5.5 c). Garcia et al. [168] investigated the defect dynamics on a sinusoidal surface using the PFC method with finite difference solver. a b d c Figure 5.5. Phase field crystal simulations on a curved surface. (a) PFC simulation of the classical Thomson problem [147]. (b) PFC simulation on the surface of the Stanford bunny [147]. (c) Defects pattern on the catenoid surface with PFC simulation [148]. (d) Defect dynamics on a sinusoidal surface using PFC [168]. 75 5.4. FEM implementation of PFC model FEM is a powerful tool to solve the PFC model on flat and curved surfaces. In this section, we will first summarize the weak forms of the governing equations of Eq. (5.9) and Eq. (5.10) and then introduce the implementation of the numerical scheme in the open source FEM software, FEniCS Project [169]. For a given function space V, the weak form for Eq. (5.9) is to search the solution  , , u V  V  V satisfying    t qdx   q  dx  0,   vdx  2 uvdx   u  vdx      3 vdx  0, (5.11)      uwdx       wdx  0, for all test functions q, v, w V . As discussed in the previous section, the higher order spatial derivatives make the stable time step for explicit time integration impractically small. Therefore, we will utilize the implicit middle point formula to deal with the time integral. The final weak form of Eq. (5.9) can be written as n1  n  qdx   q   1  n1   n dx  0,  t  2   n1vdx  2 un1vdx   un1  vdx    n1  n31 vdx  0, (5.12)       u n1wdx   n1  wdx  0,  where t  t n1  t n ,  n and  n1 denote the values at time t n and t n1 , respectively. Equation (5.12) defines a nonlinear algebra equations associated with n1 , n1 , un1 with the input n , n , un , which will be solved with the Newton method. 76 Over the past few decades, FEM has become a popular standard tool for solving PDE. A lot of effort has been devoted to developing open source FEM software package with a variety of functions and high parallel efficiency, e.g. Deal.II and FEniCS project. In the current thesis, FEniCS is used as a platform to implement the FEM scheme for the PFC method. FEniCS project is a collection of free software components for the automated, efficient solution of partial differential equations (PDE) [169], which has been widely applied to problems in fluid dynamics, surface water wave, solid mechanics and phase filed. A unique feature of FEniCS is the design and application of the Unified Form Language (UFL) [170], which is a domain-specific language, embedded in Python for specifying finite element discretizations of differential equations in terms of finite element variational forms. The UFL is relatively straightforward and easy to use, as illustrated by the following example of solving Poisson’s equation in FEniCS [170]. The continuum equation for Poisson’s equation with boundary conditions can be written as  divgrad u   f in  , u  u0 on D ,   n u  g on N , (5.13) where D  N   is a partitioning of the boundary of  into Dirichlet and Neumann boundaries. The standard H1-conforming finite element discretization of Eq. (5.13) can be expressed as: find u Vh such that au, v    grad u  grad vdx   fvdx   gvdx  Lv  , (5.14)   N for all v  Vˆh , where Vh is a continuous piecewise polynomial trial space incorporating the Dirichlet boundary conditions on D and Vˆh is a continuous piecewise polynomial test space with zero trace on D . The UFL code corresponding to the weak form defined in 77 Eq. (5.14) is shown in Fig. 5.6. It is worth noting that the derivative operator on a curved manifold has been embedded in FEniCS. Therefore, the same UFL code can work for both flat and curved surfaces, which significantly simplifies the coding effort. Figure 5.6. UFL code of FEM implementation for Poisson equation in FEniCS [170]. 5.5. Preliminary results on using PFC in curved graphene design To test our FEM code, we will first simulate the typical phase field crystal patterns shown in the phase diagram (Fig. 5.5b). A square domain with length of 19 is used in the simulations, roughly containing 3 wavelengths. The parameters in the PFC model are the same as Ref. [158]. Three typical patterns with different average densities   0.05,0.2,0.6 are simulated, corresponding to strip, triangle, and homogenous states, respectively. The initial density distribution is the flat state with a small Gaussian random perturbation, such that 0    N  ,  2  , where N  ,  2  is the Gaussian random variable with mean value  and standard derivation  and  is a small number compared with  . The final configurations for each simulation at t=1000 are summarized in Fig. 5.7, which are consistent with the theoretical predictions and previous numerical simulations. 78 a b c Figure 5.7. FEM simulations of typical PFC patterns. (a) A strip pattern for   0.05 , (b) a triangular pattern for   0.2 , and (c) a homogenous pattern for   0.6 . a b d c Figure 5.8. A flow chart for designing a curved graphene sheet via the PFC method. (a) The initial curved manifold. (b) A continuum triangular pattern on the curved surface generated by the PFC method. (c) A discrete triangular lattice network according to the triangular pattern. (d) The generated graphene atomic structure through voronoi construction from the triangular network. 79 A flow chart of employing the PFC method to design a curved 3D graphene structure is shown in Fig. 5.8, taking a sinusoidal graphene sheet as an example. The PFC simulation will be first performed on a given curved manifold (Fig. 5.8a) to search for the minimum energy triangular pattern (Fig. 5.8b). The wave crests are indentified as atoms so that the continuum triangular pattern will be transformed into a discrete triangular lattice network (Fig. 5.8c). The final graphene atomic structure (Fig. 5.8d) can then be obtained from voronoi construction based on the triangular network. For a general curved surface, including the previously discussed nanoporous graphene, the atomic structure can be generated based on a similar procedure. Once we obtain the atomic structure of the graphene sheet, we can run MD simulations to test the real configuration of the atomic structure with realistic atomic potential, such as AIREBO [46], and to investigate the mechanical and electronic properties of the obtained curved graphene. Taking the sinusoidal graphene as example, a uniaxial tensile MD simulation is first carried out to study the elastic modulus and strength of the material. Periodical boundary conditions are applied to both horizontal and vertical directions. A tensile loading with strain rate of 109/s is applied by deforming the simulation box in the horizontal direction via the command fix/deform in LAMMPS. During the simulations, the NVT ensemble is adopted to maintain the temperature at 1 K using the Nose-Hoover thermostat [56]. The stress-strain curves are depicted in Fig. 5.9a, showing a strong size dependent fracture behavior. The elastic modulus of this curved graphene is almost independent of the size and is found to be around 320 GPa by fitting the initial elastic deformation of the graphene. Typical snap shots of the stress (σxx) contours in the sinusoidal graphene sheet with different simulation size can be found in 80 Fig. 5.9b-e. The unit cell is shown together with its three periodical images. During the tensile simulation, a crack first initiates at a bond shared by a heptagon and a hexagon. For the unit cell sinusoidal graphene sheet, it is noted that catastrophic failure does not immediately follow crack nucleation (Fig. 5.9b), unlike a conventional brittle material. The reasons for this may be attributed to the small size of the sample and curved graphene geometry. On the contrary, brittle failures are found in larger sinusoidal graphene sheet (Fig. 5.9d). Another interesting phenomenon is that the sinusoidal graphene sheets remain wrinkled even after failure (Fig. 5.9c and Fig. 5.9e). a b d c e -80 80 GPa Figure 5.9. Uniaxial tensile simulation of a sinusoidal graphene sheet with molecular dynamics. (a) Stress-strain curves of the sinusoidal graphene with different unit cell size under the uniaxial tensile loading. (b)-(c) Top and front views of a snapshot of stress (σxx) contours under applied strain ε=0.15 for one unit cell with length L=4 nm. (d)-(e) Top and front views of a snapshot of stress contours under applied strain ε=0.15 for four basic unit cell with length L=16 nm. (The scale bars are 4 nm) 81 a b 8Å -8 c d e ε=0.020 ε=0.030 50 GPa 8Å ε=0.050 -50 -8 ε=0.059 82 Figure 5.10. Edge crack propagation in a nanostrip of sinusoidal and flat graphene. (a) Stress-strain curves for flat and sinusoidal graphene sheets. The dot-dashed line represents the result of a larger sinusoidal graphene strip with length of 160 nm and width of 80 nm. (b) A nanostrip of sinusoidal graphene with an edge crack, the color showing out-of plane deformation. (c) Crack tip and atomic orientation of zigzag, armchair and sinusoidal graphene. (d) A top view of stress (σyy) contours near the crack tip during the tensile deformation. (e) Contours of out-of-plane displacements from a perspective view near the crack tip during the tensile deformation. The scale bars are 10 nm and 1 nm in (b) and (d), respectively. As another illustration, we considered crack propagation in a nanostrip of sinusoidal graphene. The strip is selected to be 80 nm long and 40 nm wide, with a 20 nm long edge crack (Fig. 5.10b). Periodical boundary condition is set up in the vertical direction of the nanostrip. A tensile loading with strain rate of 109/s is applied by deforming the simulation box. During the simulations, NVT ensemble is adopted to maintain the temperature at 10 K using the Nose-Hoover thermostat [56]. From the stress strain curve, the toughness of the sinusoidal graphene is estimated to be around 32 J/m2, which is about three times that of a flat graphene (i.e. 11.4 J/m2 and 13 J/m2 for zigzag and armchair orientation, respectively). In order to exclude possible size effects on the toughness calculation, we also simulate crack propagation in a lager strip of sinusoidal graphene, with length equal to 160 nm and width equal to 80 nm, the result of which confirms the toughness value predicted using the smaller sample. More interestingly, the crack tip advances earlier in the sinusoidal graphene (around 0.023) compared to the flat graphene sheets. However, the structure remains stable so that the stress in the sample 83 keeps increasing even after the crack starts to advance. Two lattice bonds are broken at the strain of 0.03. It is noted that the crack gradually propagates until the applied strain increases to 0.059, when the heptagon ring is broken. A distinct feature of crack propagation in the curve graphene is that the local curvature keeps changing when the crack tip moves (Fig. 5.10e). It seems that this 3D curvature variation strongly influences the effective stress intensity and bond orientation near the crack tip, which may be a possible explanation for the observed toughness enhancement in sinusoidal graphene. It should be noted that the pre-stress induced by topological defects may affect both crack nucleation and propagation in the curved graphene. 5.6. Summary In this chapter, we have reviewed some previous studies related to curved graphene from both experimental and simulated points of view. It is believed that the defect and curvature effect in graphene will significantly influence the mechanical and electrical properties of graphene, such as using as a high capacity anode in Lithium ion batteries. Recent developments on the fabrication of 3D nanoporous graphene that keeps excellent electronic properties of single layer graphene have provided a promising way to design graphene based 3D electronic devices. These progresses in experimental studies are calling for more efforts in theoretical modeling and simulations on the mechanical and physical properties of the curved graphene structure to provide guidelines in real devices design. Although several attempts have been made in numerically designing curved graphene, such as Monte Carlo simulations and sphere packing mesh method, it is still challenging to generate a desired graphene structure on a given curved manifold. In this 84 regard, the phase field crystal (PFC) method is promising in searching for the minimum energy triangular pattern by tracking the dynamical diffusion of an effective continuum density function on a surface. The governing equation of the PFC model can be efficiently solved with FEM, which can be quite flexible in dealing with complicated geometries. We have implemented the FEM scheme of the PFC model in the open source software FEniCS, which is relatively easy to use and of high parallel efficiency with PETSc as the linear equation solver. Some benchmark simulations have been conducted for the phase diagram of the PFC model using the FEM code in FEniCS, which are in good agreement with previous theoretical and numerical simulations. We have employed the PFC method to create a curved graphene sheet conforming to a sinusoidal surface. The PFC simulation is first performed to search for the minimum energy triangular pattern. The wave crests are indentified as atoms so that the continuum triangular pattern is transformed into a discrete triangular lattice network. The final graphene atomic structure is obtained from voronoi construction based on the triangular network. It should be pointed out that for general curved surface, including the previously discussed nanoporous graphene, a similar procedure can be followed to generate the atomic structure. For further illustration of potential applications of this research, we further study the mechanical properties of a generated sinusoidal graphene sheet. It is demonstrated that the fracture energy of the sinusoidal graphene can be 3 times as high as that of a flat graphene sheet. 85 Chapter 6 Conclusions and Perspectives 6.1. Concluding remarks Using large-scale atomistic simulations and continuum models, we have investigated the fracture and wrinkling in graphene with topological defects. The scientific findings for each chapter are summarized as follows. In Chapter 2, flaw insensitive fracture is found in nanocrystalline graphene with a critical size around 20 nm, which is consistent with previous theoretical predictions. Our simulations reveal that micro-cracks nucleate randomly at intrinsic defects along grain boundaries and coalesce to form a big crack, eventually leading to catastrophic fracture. The key finding is that the microcrack nucleation and coalescence are not always induced by nor associated with pre-existing geometrical flaws such as a hole or notch in the system. We also extract the effective fracture energy (8 J/m2) for the nanocrystalline graphene with a 2 nm average grain size. In Chapter 3, Our MD simulations on dislocation interacting with an edge crack in a 40 nm wide graphene strip reveal that a single dislocation can increase or decrease the effective stress intensity factor by 30% depending on the orientation of the dislocation with respect to the crack. It is also interesting to note that MD simulation results lie between those predicted from the conventional linear elastic crack and the Dugdale cohesive crack model. We also carry out MD simulations for crack propagation along grain boundaries, which also have significant effect on the crack nucleation and propagation. 86 In Chapter 4, we employ large-scale atomistic simulations and continuum modeling to analyze the defects controlled wrinkling in graphene. Both atomistic and continuum simulations indicate that, due to the atomic thickness of graphene, even a single disclination/dislocation can cause significant out-of-plane wrinkling. The driving force for such wrinkling comes from the relaxation of in-plane stress and elastic energy near the defect core. Comparison with atomistic simulations indicates that the proposed model, with only three parameters (i.e., bond length, stretching modulus and bending stiffness), is capable of accurately predicting the atomic scale wrinkles near disclination/dislocation cores while also capturing the large scale graphene configurations under specific defect distributions such as those leading to a sinusoidal surface ruga or a catenoid funnel. In Chapter 5, a numerical scheme based on finite element method (FEM) is implemented in the open source software FEniCS for solving the governing equation of the phase field crystal model. We apply the phase field crystal (PFC) method to search for a triangular lattice pattern with the lowest energy on a given curved surface, which then serves as a good approximation of the graphene lattice structure conforming to that surface. The preliminary studies for the fracture of a sinusoidal graphene sheet generaed with PFC show a significant enhancement in the fracture toughness compared to the flat grpahene. 6.2. Outlook of future research It has been shown that the physical and mechanical properties of grpahene can be actively tailored through topological defects design. The studies in the current thesis have been focused on the fundamental mechanical properties of graphene with defects, with 87 the hope of providing some guidelines for real graphene based device design. Along this line of research, there are still many interesting problems, some of which are listed as follows. (1) For crack interaction with grain boundaries in graphene, only a few typical grain boundaries have been discussed in current thesis. It will be interesting to systematically simulate the crack propagation behavior and toughness for a wider range of grain boundaries. Mixed mode failure is another important phenomenon to be investigated. (2) Our preliminary studies on sinusoidal graphene Ruga have shown a promising toughening effect due to the distributed topological defects. The elastic modulus associated with curved graphene sheet is substantially decreased, which is a tradeoff for toughness enhancement. A complete map of the elastic and toughness values for different sinusoidal graphene Ruga is needed through large scale MD simulations in order to provide a balance between high modulus and toughness. (3) Curved graphene will tend to fail at bonds with high pre-stress, which is usually on a heptagon ring. For a given wavelength, although a higher wave amplitude will potentially give rise to better stretchability, it may also induce a higher pre-stress, which will accelerate the failure of the structure. Therefore, there may be an optimum aspect ratio (amplitude over wavelength) to maximize the toughness of a sinusoidal graphene Ruga. (4) Nanoporous graphene is a recent progress on designing graphene-based device, the 3D structure of which makes it easier to be manipulated. Combining phase field and phase field crystal method will generate realistic nanoporous graphene structures. It will 88 be interesting to study the mechanical and thermal properties of these structures to provide a solid foundation for real device design. (5) Although the fundamental mechanisms of the mechanical properties of graphene with defects revealed in this thesis are generally correct, we note that some values from specific simulations, like the strength and toughness of the graphene structures may vary with the interatomic potential used in the current thesis. It will be of great interest to employ more accurate potentials in future studies. 89 Bibliography [1] K. S. Novoselov, A. K. Geim, S. Morozov, D. Jiang, Y. Zhang, S. Dubonos, I. Grigorieva, A. Firsov, Electric field effect in atomically thin carbon films. science 306, 666-669 (2004). [2] K. Novoselov, A. K. Geim, S. Morozov, D. Jiang, M. K. I. Grigorieva, S. Dubonos, A. Firsov, Two-dimensional gas of massless Dirac fermions in graphene. nature 438, 197-200 (2005). [3] Y. Zhang, Y.-W. Tan, H. L. Stormer, P. Kim, Experimental observation of the quantum Hall effect and Berry's phase in graphene. Nature 438, 201-204 (2005). [4] C. Lee, X. Wei, J. W. Kysar, J. Hone, Measurement of the elastic properties and intrinsic strength of monolayer graphene. science 321, 385-388 (2008). [5] G. Eda, G. Fanchini, M. Chhowalla, Large-area ultrathin films of reduced graphene oxide as a transparent and flexible electronic material. Nature nanotechnology 3, 270-274 (2008). [6] K. S. Kim, Y. Zhao, H. Jang, S. Y. Lee, J. M. Kim, K. S. Kim, J.-H. Ahn, P. Kim, J.-Y. Choi, B. H. Hong, Large-scale pattern growth of graphene films for stretchable transparent electrodes. Nature 457, 706-710 (2009). [7] M. C. Lemme, T. J. Echtermeyer, M. Baus, H. Kurz, A graphene field-effect device. arXiv preprint cond-mat/0703208, (2007). [8] F. Xia, D. B. Farmer, Y.-m. Lin, P. Avouris, Graphene field-effect transistors with high on/off current ratio and large transport band gap at room temperature. Nano letters 10, 715-718 (2010). [9] J. S. Bunch, A. M. Van Der Zande, S. S. Verbridge, I. W. Frank, D. M. Tanenbaum, J. M. Parpia, H. G. Craighead, P. L. McEuen, Electromechanical resonators from graphene sheets. Science 315, 490-493 (2007). [10] X. Li, W. Cai, J. An, S. Kim, J. Nah, D. Yang, R. Piner, A. Velamakanni, I. Jung, E. Tutuc, Large-area synthesis of high-quality and uniform graphene films on copper foils. Science 324, 1312-1314 (2009). [11] O. V. Yazyev, S. G. Louie, Electronic transport in polycrystalline graphene. Nature materials 9, 806-809 (2010). [12] J. Lahiri, Y. Lin, P. Bozkurt, I. I. Oleynik, M. Batzill, An extended defect in graphene as a metallic wire. Nature nanotechnology 5, 326-329 (2010). [13] Q. Yu, L. A. Jauregui, W. Wu, R. Colby, J. Tian, Z. Su, H. Cao, Z. Liu, D. Pandey, D. Wei, Control and characterization of individual grains and grain boundaries in graphene grown by chemical vapour deposition. Nature materials 10, 443-449 (2011). [14] P. Y. Huang, C. S. Ruiz-Vargas, A. M. van der Zande, W. S. Whitney, M. P. Levendorf, J. W. Kevek, S. Garg, J. S. Alden, C. J. Hustedt, Y. Zhu, Grains and grain boundaries in single-layer graphene atomic patchwork quilts. Nature 469, 389-392 (2011). 90 [15] R. Grantab, V. B. Shenoy, R. S. Ruoff, Anomalous strength characteristics of tilt grain boundaries in graphene. Science 330, 946-948 (2010). [16] Y. Wei, J. Wu, H. Yin, X. Shi, R. Yang, M. Dresselhaus, The nature of strength enhancement and weakening by pentagon–heptagon defects in graphene. Nature materials 11, 759-763 (2012). [17] C. S. Ruiz-Vargas, H. L. Zhuang, P. Y. Huang, A. M. van der Zande, S. Garg, P. L. McEuen, D. A. Muller, R. G. Hennig, J. Park, Softened elastic response and unzipping in chemical vapor deposition graphene membranes. Nano letters 11, 2259-2263 (2011). [18] K. Kim, V. I. Artyukhov, W. Regan, Y. Liu, M. Crommie, B. I. Yakobson, A. Zettl, Ripping graphene: preferred directions. Nano letters 12, 293-297 (2011). [19] J. M. Carlsson, L. M. Ghiringhelli, A. Fasolino, Theory and hierarchical calculations of the structure and energetics of [0001] tilt grain boundaries in graphene. Physical Review B 84, 165423 (2011). [20] J. Zhang, J. Zhao, J. Lu, Intrinsic strength and failure behaviors of graphene grain boundaries. ACS nano 6, 2704-2711 (2012). [21] G.-H. Lee, R. C. Cooper, S. J. An, S. Lee, A. van der Zande, N. Petrone, A. G. Hammerberg, C. Lee, B. Crawford, W. Oliver, High-strength chemical-vapor– deposited graphene and grain boundaries. Science 340, 1073-1076 (2013). [22] H. I. Rasool, C. Ophus, W. S. Klug, A. Zettl, J. K. Gimzewski, Measurement of the intrinsic strength of crystalline and polycrystalline graphene. Nature communications 4, (2013). [23] A. A. Balandin, Thermal properties of graphene and nanostructured carbon materials. Nature materials 10, 569-581 (2011). [24] A. Bagri, S.-P. Kim, R. S. Ruoff, V. B. Shenoy, Thermal transport across twin grain boundaries in polycrystalline graphene from nonequilibrium molecular dynamics simulations. Nano letters 11, 3917-3921 (2011). [25] A. W. Robertson, C. S. Allen, Y. A. Wu, K. He, J. Olivier, J. Neethling, A. I. Kirkland, J. H. Warner, Spatial control of defect creation in graphene at the nanoscale. Nature communications 3, 1144 (2012). [26] M. T. Lusk, L. D. Carr, Nanoengineering defect structures on graphene. Physical review letters 100, 175503 (2008). [27] S. Zhang, S. L. Mielke, R. Khare, D. Troya, R. S. Ruoff, G. C. Schatz, T. Belytschko, Mechanics of defects in carbon nanotubes: atomistic and multiscale simulations. Physical Review B 71, 115403 (2005). [28] R. Khare, S. L. Mielke, J. T. Paci, S. Zhang, R. Ballarini, G. C. Schatz, T. Belytschko, Coupled quantum mechanical/molecular mechanical modeling of the fracture of defective carbon nanotubes and graphene sheets. Physical Review B 75, 075412 (2007). [29] S. S. Terdalkar, S. Huang, H. Yuan, J. J. Rencis, T. Zhu, S. Zhang, Nanoscale fracture in graphene. Chemical Physics Letters 494, 218-222 (2010). [30] X. Huang, H. Yang, A. C. van Duin, K. J. Hsia, S. Zhang, Chemomechanics control of tearing paths in graphene. Physical Review B 85, 195453 (2012). [31] H. Gao, B. Ji, I. L. Jäger, E. Arzt, P. Fratzl, Materials become insensitive to flaws at nanoscale: lessons from nature. Proceedings of the national Academy of Sciences 100, 5597-5600 (2003). 91 [32] O. V. Yazyev, S. G. Louie, Topological defects in graphene: Dislocations and grain boundaries. Physical Review B 81, 195420 (2010). [33] Y. Liu, B. I. Yakobson, Cones, pringles, and grain boundary landscapes in graphene topology. Nano letters 10, 2178-2183 (2010). [34] T.-H. Liu, G. Gajewski, C.-W. Pao, C.-C. Chang, Structure, energy, and structural transformations of graphene grain boundaries from atomistic simulations. Carbon 49, 2306-2317 (2011). [35] C.-W. Pao, T.-H. Liu, C.-C. Chang, D. J. Srolovitz, Graphene defect polarity dynamics. Carbon 50, 2870-2876 (2012). [36] O. Lehtinen, S. Kurasch, A. Krasheninnikov, U. Kaiser, Atomic scale study of the life cycle of a dislocation in graphene from birth to annihilation. Nature communications 4, (2013). [37] J. H. Warner, Y. Fan, A. W. Robertson, K. He, E. Yoon, G. D. Lee, Rippling graphene at the nanoscale through dislocation addition. Nano letters 13, 4937- 4944 (2013). [38] N. Levy, S. Burke, K. Meaker, M. Panlasigui, A. Zettl, F. Guinea, A. C. Neto, M. Crommie, Strain-induced pseudo–magnetic fields greater than 300 tesla in graphene nanobubbles. Science 329, 544-547 (2010). [39] V. M. Pereira, A. C. Neto, H. Liang, L. Mahadevan, Geometry, mechanics, and electronics of singular structures and wrinkles in graphene. Physical review letters 105, 156603 (2010). [40] N. N. Klimov, S. Jung, S. Zhu, T. Li, C. A. Wright, S. D. Solares, D. B. Newell, N. B. Zhitenev, J. A. Stroscio, Electromechanical properties of graphene drumheads. Science 336, 1557-1561 (2012). [41] J. S. Choi, J.-S. Kim, I.-S. Byun, D. H. Lee, M. J. Lee, B. H. Park, C. Lee, D. Yoon, H. Cheong, K. H. Lee, Friction anisotropy–driven domain imaging on exfoliated monolayer graphene. Science 333, 607-610 (2011). [42] J. Zang, S. Ryu, N. Pugno, Q. Wang, Q. Tu, M. J. Buehler, X. Zhao, Multifunctionality and control of the crumpling and unfolding of large-area graphene. Nature materials 12, 321-325 (2013). [43] M. P. Allen, D. J. Tildesley, Computer simulation of liquids. (1987). [44] J. Tersoff, New empirical approach for the structure and energy of covalent systems. Physical Review B 37, 6991 (1988). [45] D. W. Brenner, O. A. Shenderova, J. A. Harrison, S. J. Stuart, B. Ni, S. B. Sinnott, A second-generation reactive empirical bond order (REBO) potential energy expression for hydrocarbons. Journal of Physics: Condensed Matter 14, 783 (2002). [46] S. J. Stuart, A. B. Tutein, J. A. Harrison, A reactive potential for hydrocarbons with intermolecular interactions. The Journal of chemical physics 112, 6472-6486 (2000). [47] N. Marks, Generalizing the environment-dependent interaction potential for carbon. Physical Review B 63, 035401 (2000). [48] J. H. Los, L. M. Ghiringhelli, E. J. Meijer, A. Fasolino, Improved long-range reactive bond-order potential for carbon. I. Construction. Physical Review B 72, 214102 (2005). 92 [49] A. C. Van Duin, S. Dasgupta, F. Lorant, W. A. Goddard, ReaxFF: a reactive force field for hydrocarbons. The Journal of Physical Chemistry A 105, 9396-9409 (2001). [50] T. Liang, T.-R. Shan, Y.-T. Cheng, B. D. Devine, M. Noordhoek, Y. Li, Z. Lu, S. R. Phillpot, S. B. Sinnott, Classical atomistic simulations of surfaces and heterogeneous interfaces with the charge-optimized many body (COMB) potentials. Materials Science and Engineering: R: Reports 74, 255-279 (2013). [51] R. Perriot, X. Gu, Y. Lin, V. V. Zhakhovsky, I. I. Oleynik, Screened environment-dependent reactive empirical bond-order potential for atomistic simulations of carbon materials. Physical Review B 88, 064101 (2013). [52] R. Z. Khaliullin, H. Eshet, T. D. Kühne, J. Behler, M. Parrinello, Nucleation mechanism for the direct graphite-to-diamond phase transition. Nature materials 10, 693-697 (2011). [53] S. Plimpton, Fast parallel algorithms for short-range molecular dynamics. Journal of computational physics 117, 1-19 (1995). [54] O. Shenderova, D. Brenner, A. Omeltchenko, X. Su, L. Yang, Atomistic modeling of the fracture of polycrystalline diamond. Physical Review B 61, 3877 (2000). [55] H. J. Berendsen, J. P. M. Postma, W. F. van Gunsteren, A. DiNola, J. Haak, Molecular dynamics with coupling to an external bath. The Journal of chemical physics 81, 3684-3690 (1984). [56] W. G. Hoover, Canonical dynamics: equilibrium phase-space distributions. Physical Review A 31, 1695 (1985). [57] S. Nosé, A unified formulation of the constant temperature molecular dynamics methods. The Journal of Chemical Physics 81, 511-519 (1984). [58] T. Schneider, E. Stoll, Dynamics of the Sine-Gordon Chain. Physical Review Letters 41, 1429 (1978). [59] H. C. Andersen, Molecular dynamics simulations at constant pressure and/or temperature. The Journal of chemical physics 72, 2384-2393 (1980). [60] M. Parrinello, A. Rahman, Polymorphic transitions in single crystals: A new molecular dynamics method. Journal of Applied physics 52, 7182-7190 (1981). [61] M. Diab, T. Zhang, R. Zhao, H. Gao, K.-S. Kim, Ruga mechanics of creasing: from instantaneous to setback creases. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Science 469, 20120753 (2013). [62] J. Bai, X. Zhong, S. Jiang, Y. Huang, X. Duan, Graphene nanomesh. Nature nanotechnology 5, 190-194 (2010). [63] W. Oswald, Z. Wu, Energy gaps in graphene nanomeshes. Physical Review B 85, 115431 (2012). [64] D. Cohen-Tanugi, J. C. Grossman, Water desalination across nanoporous graphene. Nano letters 12, 3602-3608 (2012). [65] G. F. Schneider, S. W. Kowalczyk, V. E. Calado, G. Pandraud, H. W. Zandbergen, L. M. Vandersypen, C. Dekker, DNA translocation through graphene nanopores. Nano letters 10, 3163-3167 (2010). [66] H. W. C. Postma, Rapid sequencing of individual DNA molecules in graphene nanogaps. Nano letters 10, 420-425 (2010). 93 [67] Z. P. Bazant, J. Planas, Fracture and size effect in concrete and other quasibrittle materials. (CRC press, 1997), vol. 16. [68] S. J. Bennison, B. R. Lawn, Flaw tolerance in ceramics with rising crack resistance characteristics. Journal of materials science 24, 3169-3175 (1989). [69] S. Kumar, M. Haque, H. Gao, Notch insensitive fracture in nanoscale thin films. Applied Physics Letters 94, 253104 (2009). [70] S. Kumar, X. Li, A. Haque, H. Gao, Is stress concentration relevant for nanocrystalline metals? Nano letters 11, 2510-2516 (2011). [71] Q. Lu, N. Marks, G. C. Schatz, T. Belytschko, Nanoscale fracture of tetrahedral amorphous carbon by molecular dynamics: Flaw size insensitivity. Physical Review B 77, 014109 (2008). [72] Z. Qin, M. J. Buehler, Flaw tolerance of nuclear intermediate filament lamina under extreme mechanical deformation. Acs Nano 5, 3034-3042 (2011). [73] T. Giesa, M. Arslan, N. M. Pugno, M. J. Buehler, Nanoconfinement of spider silk fibrils begets superior strength, extensibility, and toughness. Nano letters 11, 5038-5046 (2011). [74] S. Keten, Z. Xu, B. Ihle, M. J. Buehler, Nanoconfinement controls stiffness, strength and mechanical toughness of [beta]-sheet crystals in silk. Nature materials 9, 359-367 (2010). [75] Z.-D. Sha, Q.-X. Pei, V. Sorkin, P. S. Branicio, Y.-W. Zhang, H. Gao, On the notch sensitivity of CuZr metallic glasses. Applied Physics Letters 103, 081903 (2013). [76] Z. Sha, L. He, Q. Pei, H. Pan, Z. Liu, Y. Zhang, T. Wang, On the notch sensitivity of CuZr nanoglass. Journal of Applied Physics 115, 163507 (2014). [77] X. W. Gu, Z. Wu, Y.-W. Zhang, D. J. Srolovitz, J. R. Greer, Microstructure versus Flaw: Mechanisms of Failure and Strength in Nanostructures. Nano letters 13, 5703-5709 (2013). [78] H. Gao, S. Chen, Flaw tolerance in a thin strip under tension. Journal of applied mechanics 72, 732-737 (2005). [79] J. C. Meyer, A. K. Geim, M. Katsnelson, K. Novoselov, T. Booth, S. Roth, The structure of suspended graphene sheets. Nature 446, 60-63 (2007). [80] A. Fasolino, J. Los, M. I. Katsnelson, Intrinsic ripples in graphene. Nature materials 6, 858-861 (2007). [81] M. J. Buehler, H. Yao, H. Gao, B. Ji, Cracking and adhesion at small scales: atomistic and continuum studies of flaw tolerant nanostructures. Modelling and Simulation in Materials Science and Engineering 14, 799 (2006). [82] H. Zhao, N. Aluru, Temperature and strain-rate dependent fracture strength of graphene. Journal of Applied Physics 108, 064321 (2010). [83] T. L. Anderson, Fracture mechanics: fundamentals and applications. (CRC press, 2005). [84] Y. Jin, F. Yuan, Nanoscopic modeling of fracture of 2D graphene systems. Journal of nanoscience and nanotechnology 5, 601-608 (2005). [85] J. P. Hirth, J. Lothe, Theory of dislocations. (1982). [86] P. Zhang, L. Ma, F. Fan, Z. Zeng, C. Peng, P. E. Loya, Z. Liu, Y. Gong, J. Zhang, X. Zhang, Fracture toughness of graphene. Nature communications 5, (2014). 94 [87] J. Li, C.-Z. Wang, J.-P. Chang, W. Cai, V. V. Bulatov, K.-M. Ho, S. Yip, Core energy and Peierls stress of a screw dislocation in bcc molybdenum: A periodic- cell tight-binding study. Physical Review B 70, 104113 (2004). [88] E. Ertekin, D. Chrzan, M. S. Daw, Topological description of the Stone-Wales defect formation energy in carbon nanotubes and graphene. Physical Review B 79, 155421 (2009). [89] I.-H. Lin, R. Thomson, Cleavage, dislocation emission, and shielding for cracks under general loading. Acta Metallurgica 34, 187-206 (1986). [90] T. K. Bhandakkar, A. C. Chng, W. Curtin, H. Gao, Dislocation shielding of a cohesive crack. Journal of the Mechanics and Physics of Solids 58, 530-541 (2010). [91] J. Song, W. Curtin, T. Bhandakkar, H. Gao, Dislocation shielding and crack tip decohesion at the atomic scale. Acta Materialia 58, 5933-5940 (2010). [92] D. Dugdale, Yielding of steel sheets containing slits. Journal of the Mechanics and Physics of Solids 8, 100-104 (1960). [93] J. Li, Disclination model of high angle grain boundaries. Surface Science 31, 12- 26 (1972). [94] Y. Cheng, Z.-H. Jin, Y. Zhang, H. Gao, On intrinsic brittleness and ductility of intergranular fracture along symmetrical tilt grain boundaries in copper. Acta Materialia 58, 2293-2299 (2010). [95] J. R. Rice, Dislocation nucleation from a crack tip: an analysis based on the Peierls concept. Journal of the Mechanics and Physics of Solids 40, 239-271 (1992). [96] J. J. Möller, E. Bitzek, Fracture toughness and bond trapping of grain boundary cracks. Acta Materialia 73, 1-11 (2014). [97] C.-K. Wang, H. B. Chew, K.-S. Kim, in MRS Proceedings. (Cambridge Univ Press, 2011), vol. 1297, pp. mrsf10-1297-p1202-1201. [98] G. Xu, M. Demkowicz, Healing of Nanocracks by Disclinations. Physical review letters 111, 145501 (2013). [99] D. Nelson, L. Peliti, Fluctuations in membranes with crystalline and hexatic order. Journal de physique 48, 1085-1092 (1987). [100] H. Seung, D. R. Nelson, Defects in flexible membranes with crystalline order. Physical Review A 38, 1005 (1988). [101] L. M. Zubov, Nonlinear theory of dislocations and disclinations in elastic bodies. (Springer, 1997), vol. 47. [102] L. Zubov, in Doklady Physics. (Springer, 2007), vol. 52, pp. 67-70. [103] L. Zubov, The linear theory of dislocations and disclinations in elastic shells. Journal of Applied Mathematics and Mechanics 74, 663-672 (2010). [104] Y. Huang, J. Wu, K. Hwang, Thickness of graphene and single-wall carbon nanotubes. Physical review B 74, 245413 (2006). [105] J. Wu, K. Hwang, Y. Huang, An atomistic-based finite-deformation shell theory for single-wall carbon nanotubes. Journal of the Mechanics and Physics of Solids 56, 279-292 (2008). [106] Q. Lu, M. Arroyo, R. Huang, Elastic bending modulus of monolayer graphene. Journal of Physics D: Applied Physics 42, 102002 (2009). 95 [107] X. Shi, B. Peng, N. M. Pugno, H. Gao, Stretch-induced softening of bending rigidity in graphene. Applied Physics Letters 100, 191913 (2012). [108] J. Wu, Y. Wei, Grain misorientation and grain-boundary rotation dependent mechanical properties in polycrystalline graphene. Journal of the Mechanics and Physics of Solids 61, 1421-1432 (2013). [109] S. Chen, D. Chrzan, Continuum theory of dislocations and buckling in graphene. Physical Review B 84, 214103 (2011). [110] T. Mura, Micromechanics of defects in solids. (Springer, 1987), vol. 3. [111] J. Dervaux, M. B. Amar, Morphogenesis of growing soft tissues. Physical review letters 101, 068101 (2008). [112] H. Liang, L. Mahadevan, The shape of a long leaf. Proceedings of the National Academy of Sciences 106, 22049-22054 (2009). [113] B. Li, Y.-P. Cao, X.-Q. Feng, H. Gao, Surface wrinkling of mucosa induced by volumetric growth: theory, simulation and experiment. Journal of the Mechanics and Physics of Solids 59, 758-774 (2011). [114] B. Li, F. Jia, Y.-P. Cao, X.-Q. Feng, H. Gao, Surface wrinkling patterns on a core- shell soft sphere. Physical review letters 106, 234301 (2011). [115] Y. Klein, E. Efrati, E. Sharon, Shaping of elastic sheets by prescription of non- Euclidean metrics. Science 315, 1116-1120 (2007). [116] J. Kim, J. A. Hanna, M. Byun, C. D. Santangelo, R. C. Hayward, Designing responsive buckled surfaces by halftone gel lithography. Science 335, 1201-1205 (2012). [117] E. Efrati, E. Sharon, R. Kupferman, Elastic theory of unconstrained non- Euclidean plates. Journal of the Mechanics and Physics of Solids 57, 762-775 (2009). [118] J. Dervaux, P. Ciarletta, M. Ben Amar, Morphogenesis of thin hyperelastic plates: a constitutive theory of biological growth in the Föppl–von Kármán limit. Journal of the Mechanics and Physics of Solids 57, 458-471 (2009). [119] M. Lewicka, L. Mahadevan, M. R. Pakzad, The Föppl-von Kármán equations for plates with incompatible strains. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Science 467, 402-426 (2011). [120] B. Li, Y.-P. Cao, X.-Q. Feng, H. Gao, Mechanics of morphological instabilities and surface wrinkling in soft materials: a review. Soft Matter 8, 5728-5745 (2012). [121] J. D. Eshelby, The determination of the elastic field of an ellipsoidal inclusion, and related problems. Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences 241, 376-396 (1957). [122] Y. Wei, B. Wang, J. Wu, R. Yang, M. L. Dunn, Bending rigidity and Gaussian bending stiffness of single-layered graphene. Nano letters 13, 26-30 (2012). [123] X. Zhang, K. Jiao, P. Sharma, B. Yakobson, An atomistic and non-classical continuum field theoretic perspective of elastic interactions between defects (force dipoles) of various symmetries and application to graphene. Journal of the Mechanics and Physics of Solids 54, 2304-2329 (2006). [124] A. Carpio, L. L. Bonilla, F. de Juan, M. A. Vozmediano, Dislocations in graphene. New Journal of Physics 10, 053021 (2008). 96 [125] M. Ariza, M. Ortiz, R. Serrano, Long-term dynamic stability of discrete dislocations in graphene at finite temperature. International journal of fracture 166, 215-223 (2010). [126] E. Cerda, L. Mahadevan, Geometry and physics of wrinkling. Physical review letters 90, 074302 (2003). [127] F. Liu, P. Ming, J. Li, Ab initio calculation of ideal strength and phonon instability of graphene under tension. Physical Review B 76, 064120 (2007). [128] A. Bausch, M. Bowick, A. Cacciuto, A. Dinsmore, M. Hsu, D. Nelson, M. Nikolaides, A. Travesset, D. Weitz, Grain boundary scars and spherical crystallography. Science 299, 1716-1718 (2003). [129] V. Vitelli, J. B. Lucks, D. R. Nelson, Crystallography on curved surfaces. Proceedings of the National Academy of Sciences 103, 12323-12328 (2006). [130] A. Hexemer, V. Vitelli, E. J. Kramer, G. H. Fredrickson, Monte Carlo study of crystalline order and defects on weakly curved surfaces. Physical Review E 76, 051604 (2007). [131] M. J. Bowick, L. Giomi, Two-dimensional matter: order, curvature and defects. Advances in Physics 58, 449-563 (2009). [132] W. T. Irvine, V. Vitelli, P. M. Chaikin, Pleats in crystals on curved surfaces. Nature 468, 947-951 (2010). [133] W. T. Irvine, M. J. Bowick, P. M. Chaikin, Fractionalization of interstitials in curved colloidal crystals. Nature materials 11, 948-951 (2012). [134] H. Kusumaatmaja, D. J. Wales, Defect Motifs for Constant Mean Curvature Surfaces. Physical review letters 110, 165502 (2013). [135] J. Lidmar, L. Mirny, D. R. Nelson, Virus shapes and buckling transitions in spherical shells. Physical Review E 68, 051910 (2003). [136] C. M. Funkhouser, R. Sknepnek, M. O. de la Cruz, Topological defects in the buckling of elastic membranes. Soft Matter 9, 60-68 (2013). [137] E. H. Yong, D. R. Nelson, L. Mahadevan, Elastic Platonic Shells. Physical review letters 111, 177801 (2013). [138] C. D. Modes, K. Bhattacharya, M. Warner, Disclination-mediated thermo-optical response in nematic glass sheets. Physical Review E 81, 060701 (2010). [139] C. Modes, K. Bhattacharya, M. Warner, Gaussian curvature from flat elastica sheets. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Science 467, 1121-1140 (2011). [140] K. Kawasumi, Q. Zhang, Y. Segawa, L. T. Scott, K. Itami, A grossly warped nanographene and the consequences of multiple odd-membered-ring defects. Nature chemistry, (2013). [141] C. Chuang, Y.-C. Fan, B.-Y. Jin, Generalized classification scheme of toroidal and helical carbon nanotubes. Journal of chemical information and modeling 49, 361-368 (2009). [142] C. Chuang, B.-Y. Jin, Systematics of high-genus fullerenes. Journal of chemical information and modeling 49, 1664-1668 (2009). [143] E. Biyikli, J. Liu, X. Yang, A. C. To, A fast method for generating atomistic models of arbitrarily-shaped carbon graphitic nanostructures. RSC Advances 3, 1359-1362 (2013). 97 [144] T. C. Petersen, I. K. Snook, I. Yarovsky, D. G. McCulloch, Monte Carlo based modeling of carbon nanostructured surfaces. Physical Review B 72, 125417 (2005). [145] K. Elder, M. Katakowski, M. Haataja, M. Grant, Modeling elasticity in crystal growth. Physical review letters 88, 245701 (2002). [146] R. Backofen, A. Voigt, T. Witkowski, Particles on curved surfaces: A dynamic approach by a phase-field-crystal model. Physical Review E 81, 025701 (2010). [147] R. Backofen, M. Gräf, D. Potts, S. Praetorius, A. Voigt, T. Witkowski, A Continuous Approach to Discrete Ordering on S^2. Multiscale Modeling & Simulation 9, 314-334 (2011). [148] V. Schmid, A. Voigt, Crystalline order and topological charges on capillary bridges. Soft matter, (2014). [149] H. Terrones, M. Terrones, Curved nanostructured materials. New Journal of Physics 5, 126 (2003). [150] H. Terrones, A. Mackay, The geometry of hypothetical curved graphite structures. Carbon 30, 1251-1260 (1992). [151] T. Lenosky, X. Gonze, M. Teter, V. Elser, Energetics of negatively curved graphitic carbon. (1992). [152] T. Ichihashi, Y. Ando, Pentagons, heptagons and negative curvature in graphite microtubule growth. Nature 356, 776-778 (1992). [153] Y. Zhu, S. Murali, M. D. Stoller, K. Ganesh, W. Cai, P. J. Ferreira, A. Pirkle, R. M. Wallace, K. A. Cychosz, M. Thommes, Carbon-based supercapacitors produced by activation of graphene. Science 332, 1537-1541 (2011). [154] D. Odkhuu, D. H. Jung, H. Lee, S. S. Han, S.-H. Choi, R. S. Ruoff, N. Park, Negatively curved carbon as the anode for lithium ion batteries. Carbon 66, 39-47 (2014). [155] Y. Ito, Y. Tanabe, H. J. Qiu, K. Sugawara, S. Heguri, N. H. Tu, K. K. Huynh, T. Fujita, T. Takahashi, K. Tanigaki, High‐Quality Three‐Dimensional Nanoporous Graphene. Angewandte Chemie International Edition 53, 4822-4826 (2014). [156] H. Emmerich, H. Löwen, R. Wittkowski, T. Gruhn, G. I. Tóth, G. Tegze, L. Gránásy, Phase-field-crystal models for condensed matter dynamics on atomic length and diffusive time scales: an overview. Advances in Physics 61, 665-743 (2012). [157] K. Elder, M. Grant, Modeling elastic and plastic deformations in nonequilibrium processing using phase field crystals. Physical Review E 70, 051605 (2004). [158] P. Galenko, H. Gomez, N. Kropotin, K. Elder, Unconditionally stable method and numerical solution of the hyperbolic phase-field crystal equation. Physical Review E 88, 013310 (2013). [159] J. Mellenthin, A. Karma, M. Plapp, Phase-field crystal study of grain-boundary premelting. Physical Review B 78, 184110 (2008). [160] L. Gránásy, T. Pusztai, J. A. Warren, Modelling polycrystalline solidification using phase field theory. Journal of physics: Condensed matter 16, R1205 (2004). [161] K.-A. Wu, P. W. Voorhees, Phase field crystal simulations of nanocrystalline grain growth in two dimensions. Acta Materialia 60, 407-419 (2012). 98 [162] K. Elder, N. Provatas, J. Berry, P. Stefanovic, M. Grant, Phase-field crystal modeling and classical density functional theory of freezing. Physical Review B 75, 064107 (2007). [163] I. Singer-Loginova, H. Singer, The phase field technique for modeling multiphase materials. Reports on progress in physics 71, 106501 (2008). [164] K.-A. Wu, A. Karma, Phase-field crystal modeling of equilibrium bcc-liquid interfaces. Physical Review B 76, 184107 (2007). [165] S. Mkhonta, K. Elder, Z.-F. Huang, Exploring the complex world of two- dimensional ordering with three modes. Physical review letters 111, 035501 (2013). [166] J. Berry, N. Provatas, J. Rottler, C. W. Sinclair, Phase field crystal modeling as a unified atomistic approach to defect dynamics. Physical Review B 89, 214117 (2014). [167] R. Backofen, A. Rätz, A. Voigt, Nucleation and growth by a phase field crystal (PFC) model. Philosophical Magazine Letters 87, 813-820 (2007). [168] N. A. García, R. A. Register, D. A. Vega, L. R. Gómez, Crystallization dynamics on curved surfaces. Physical Review E 88, 012306 (2013). [169] A. Logg, G. N. Wells, DOLFIN: Automated finite element computing. ACM Transactions on Mathematical Software (TOMS) 37, 20 (2010). [170] M. S. Alnaes, A. Logg, K. B. Ølgaard, M. E. Rognes, G. N. Wells, Unified form language: a domain-specific language for weak formulations of partial differential equations. ACM Transactions on Mathematical Software (TOMS) 40, 9 (2014). 99