UFDC Home  myUFDC Home  Help 



Full Text  
PAGE 1 1 COARSE GRAINING ATOMISTIC DYNAMICS OF FRACTURE BY FINITE ELEMENT METHOD: FORM UL ATION, PARALLELIZATION AND APPLICATIONS By QIAN DENG A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFIL LMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 2011 PAGE 2 2 2011 Qian Deng PAGE 3 3 To my wife and mother PAGE 4 4 ACKNOWLEDGMENTS Though only my name appears on the cover of this dissertation, a great many people have cont ributed to its production. I owe my gratitude to all those people who have made this dissertation possible and because of whom my graduate experience has been one that I will cherish forever. My deepest gratitude is to my advisor, Dr. Youping Chen. In the past four years, she has always been there to support and give guidance to me. S he taught me how to do research and showed me the frontier of my research field Most importantly, she taught me the way to explore the unknown by myself. I have learned a grea t deal from her and I will never forget the valuable lessons she taught me. I hope that one day I would become a scholar as good as Dr. Chen. I would also like to thank Dr James Le e for his excellent guidance, caring and patience when I was in t he George Washington University I would like to thank Dr. Loc Vu Quoc for his valuable discussions with me that helped me understand my research area better I would like to thank Dr. Bhavan i V. Sankar for giving me kind help and advi ce during the past years I would also like to thank Dr. Beverly A. Sanders for helping me to develop my background in parallel computing. Special thank s go to Dr. Peter G. Ifju who was willing to participate in my final defense committee in the last semester. I would also like to thank my wife and my parents for their encouragement and unwavering support throughout my life. Their love is always providing me ins piration and is my driving force. I could not have come this far without their assistance and I want to express my deepest appreciation to them. I would also like to thank Dr. Liming Xiong, my lab mate, for his discussion s and PAGE 5 5 cooperation i n the developmen t of the computational codes in this work. I would also like to thank all of the other lab mates, Ms. Ning Zhang, Mr. Shi Li, Mr. Shengfeng Yang, Mr. Yongtao Li and Mr. Yu Hong, for providing such a thought provoking and cooperative atmosphere in our lab My special thanks also go to Mr. Eddie McCumiskey my lab mate, for his kind help and valuble suggestions in the past few yea rs. Finally, the generous financial support from NSF (National Science Foundation) DARPA ( Defense Advanced Research Projects Agency ) and DOE (Department of Energy) on this project is highly appreciated. PAGE 6 6 TABLE OF CONTENTS p age ACKNOWLEDGMENTS ................................ ................................ ................................ .. 4 LIST OF TABLES ................................ ................................ ................................ ............ 8 LIST OF FIG URES ................................ ................................ ................................ .......... 9 ABSTRACT ................................ ................................ ................................ .................... 11 CHAPTER 1 MOTIVATIONS ................................ ................................ ................................ ......... 13 2 BACKGROUND ................................ ................................ ................................ ....... 18 State of the art of numerical methods for fracture ................................ ................... 18 Continuum Methods ................................ ................................ ......................... 18 Atomistic Simulation Methods ................................ ................................ .......... 22 Multiscale Simulations Methods ................................ ................................ ....... 24 Sequential multiscale simulations: ................................ ............................. 24 Concurrent multiscale simulations: ................................ ............................ 25 Summary ................................ ................................ ................................ .......... 27 Introduction to parallel computi ng ................................ ................................ ........... 28 What is Parallel Computing? ................................ ................................ ............ 28 Why Do We Use Parallel Computing Here? ................................ ..................... 28 High Performance Computing Resources ................................ ........................ 29 Message Passing Interface (MPI) ................................ ................................ .... 30 Running a Parallel Code on High Perfor mance Computers ............................. 33 Summary ................................ ................................ ................................ .......... 33 3 A COARSE GRAINED METHODOLOGY BASED ON THE ATOMISTIC FIELD THEORY (AFT) ................................ ................................ ................................ ........ 35 Introduction ................................ ................................ ................................ ............. 35 A brief review of coarse grained (CG) models ................................ ........................ 35 An introduction to AFT ................................ ................................ ............................ 38 Decomposition of the Atomic Motion ................................ ................................ 38 The Local Density and Its Time Derivative ................................ ....................... 39 Balance Equations ................................ ................................ ........................... 41 Difference between Particles Dynamics and Continuous Field Theories .......... 44 The finite el ement implementation of AFT ................................ ............................... 45 Weak Form of the Governing Equation ................................ ............................ 45 Numerical Evaluation of Each Term in Weak Form Equation ........................... 47 Mesh Strategies ................................ ................................ ................................ 51 2D triangular element ................................ ................................ ................. 52 3D hybrid element ................................ ................................ ...................... 53 PAGE 7 7 Summary ................................ ................................ ................................ ................ 54 4 PARALLEL IMPLEMENTATION OF CG AFT METHOD ................................ .......... 59 Intro duction ................................ ................................ ................................ ............. 59 Four design spaces for the parallel code ................................ ................................ 59 The structure of the parallel code ................................ ................................ ........... 62 Domain Decomposition, GETMAP ................................ ................................ ... 63 Construct Neighborlist, ELIST ................................ ................................ .......... 63 Communication Between Processors, EXCHANGE_U ................................ .... 64 Calculation of Nodal Forces, CFORCE ................................ ............................ 65 Time Integrations, CENTRAL ................................ ................................ ........... 65 Performance of the parallel code ................................ ................................ ............ 65 Stability ................................ ................................ ................................ ............. 65 Scalability ................................ ................................ ................................ ......... 66 Summary ................................ ................................ ................................ ................ 69 5 SIMULATIONS OF 2D DYNAMIC BRITTLE FRACTURE ................................ ........ 76 Introduction ................................ ................................ ................................ ............. 76 The computational model ................................ ................................ ........................ 76 Simulation details and results ................................ ................................ ................. 78 Comparison of CG And MD Simulations ................................ ........................... 78 Large Scale FE Simulations ................................ ................................ ............. 81 Summary ................................ ................................ ................................ ................ 83 6 SIMULATIONS OF 3D D YNAMIC BRITTLE FRACTURE ................................ ........ 95 Introduction ................................ ................................ ................................ ............. 95 Comparison of CG and MD simulations ................................ ................................ .. 96 Application of the CG method in 3D brittle fracture ................................ ................. 98 3D Crack Branching Under Mixed Mode Loading ................................ ............ 98 Impact F racture of a Plate ................................ ................................ ................ 99 3D Crack Surface Evolution ................................ ................................ ........... 101 Summary ................................ ................................ ................................ .............. 102 7 SUMMARY AND FUTURE WORK ................................ ................................ .......... 113 Summary and conclusion ................................ ................................ ....................... 113 Future work ................................ ................................ ................................ ............ 115 LIST OF REFERENCES ................................ ................................ .............................. 117 BIOGRAPHICAL SKETCH ................................ ................................ .......................... 123 PAGE 8 8 LIST OF TABLES Table page 4 1 Number of elements and processors used in ST and FT ................................ .... 70 4 2 FT results (unit: second) ................................ ................................ ..................... 70 4 3 ST results (unit: second) ................................ ................................ ..................... 70 5 1 Geometric parameters and loading rates for the six models ............................... 85 5 2 K IC for models A~E ................................ ................................ ............................. 85 5 3 Time interval during which crack branching occurs and time of arrival of reflected stress waves ................................ ................................ ................................ ....... 85 PAGE 9 9 LIST OF FIGURES Figure page 1 1 Crack oscillation and branching under mode I loading ................................ ....... 16 1 2 Helical crack front under mixed mode (I+III) loading ................................ .......... 16 1 3 Fracture at different length scales ................................ ................................ ...... 16 1 4 Fracture of some single crystals ................................ ................................ ......... 17 2 1 Cohesive element and bulk element ................................ ................................ ... 34 3 1 Relationship between AFT and continuum mechanics ................................ ....... 56 3 2 Decomposition of atomic positions ................................ ................................ ..... 56 3 3 Schematic scratch for finite elements in the FE implementation of AFT ............. 56 3 4 Nodal integration and Gauss integration ................................ ............................ 57 3 5 Two options of the mesh strategy ................................ ................................ ....... 57 3 6 The triangular element which shows the weakest planes ................................ ... 57 3 7 T he primitive unit cell and slip planes of FCC crystal ................................ ......... 58 3 8 Sketch of one hybrid element ................................ ................................ ............. 58 4 1 Structure of the parallel code ................................ ................................ .............. 71 4 2 Procedures for domain decomposition ................................ ............................... 71 4 3 Interfaces between domains ................................ ................................ ............... 71 4 4 Communication between neighboring domains ................................ .................. 72 4 5 Central difference method ................................ ................................ .................. 72 4 6 Setu p of the model for the force perturbation test ................................ ............... 73 4 7 Vibration of a FCC bar predicted by the two different methods .......................... 73 4 8 Vibra tion of the FCC bar predicted with varies number of processors ................ 74 4 9 Speedup curves for FT and ST ................................ ................................ ........... 74 4 10 Efficiency curves for FT and ST ................................ ................................ .......... 75 PAGE 10 10 5 1 Simulation Model and Boundary Conditions ................................ ....................... 86 5 2 The finite element mesh of the computer model ................................ ................. 86 5 3 Comparison between FE and MD simulations results ................................ ........ 87 5 4 Average stress vs. strain curves ................................ ................................ ......... 88 5 5 Comparison of FE and MD simulations results ................................ ................... 88 5 6 Stress distribution in front of crack tip when crack start to propagate ................. 90 5 7 Stress distribution in front of crack tip when crack start to propagate ................. 90 5 8 FE resul ts of dynamic crack propagation ................................ ............................ 91 5 9 Crack speed measured for model E during the period of crack initiation and first branching ................................ ................................ ................................ ............ 93 5 10 Stress waves at the time of crack branching in models A E .............................. 93 6 1 Notched film for the first test case ................................ ................................ .... 104 6 2 Comparison of CG and MD results for the first case ................................ ......... 104 6 3 Comparison of the average stress time curves ................................ .............. 105 6 4 Computation model for 3D crack branching ................................ ...................... 105 6 5 MD simulation of crack branching under mixed mode loading ......................... 106 6 6 CG simulation of crack branching under mixed mode loading. ........................ 106 6 7 Computation model for 3D impact test ................................ .............................. 107 6 8 Fracture of the plate under the impact of a rigid ball ................................ ......... 108 6 9 Impact test for a smaller plate by MD simulation ................................ .............. 108 6 10 Computation model for the bulk material under Mode I loading ........................ 109 6 11 Crack surface evolution in bulk material ................................ ........................... 109 PAGE 11 11 Abstract of Dis sertation Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy COARSE GRAINING ATOMISTIC DYNAMICS OF FRACTURE BY FINITE ELEMENT METHOD: FORMULATION, PARALLELIZATION AND APPLICATIONS By Qian Deng December 2011 Chair: Youping Chen Major: Mechanical Engineering The accurate prediction of material behaviors is one o f the most important fields in both science and engineering communities. Understanding the mechanisms of material behaviors sometimes needs us to study the material over a wide span of length scales. In this work we present a new metho dology which is able to coarse g rain the atomistic dynamics of fracture by finite element method First bas ed on the Atomistic Field Theory (AFT) (Chen and Lee 2005, Chen et al. 2006, Chen 2009) a finite elements method with built in atomistic information is presented Then a high efficiency parallel code for large scale computation is described. This code was written in FORTRAN language and uses the standard parallel programming environment me ssage passing interface (MPI ). The pe rformance of the parallel code was tested on the supercomputer Trestle of SDSC ( San Diego Supercomputer Center ). Through the comparis on of the coarse grained (CG) simulation results with the molecular dynamics (MD) simulation results it is found that the new CG method is able to predict the crack tip stress, dynamic crack propagation and even crack branching, with results similar to th at of the atomic level molecular dynamics PAGE 12 12 simu lations Finally, both 2D (2 dimensional) and 3D (3 dimensional) dynamic fracture problems were computed through the CG method In 2D dynamic fracture simulations, the relationship between stress waves and crac k propagation s wa s studied. It is found that the stress waves reflected back from the boundary can trigger the dynamic crack branching. In 3D simulations the dynamic fractures under different loading were simulated T he largest 3D model i s composed of ove r 0.1 million elements which are equivalent to over 0.1 billion atom s. To show the performance of the parallel code in dealing with large number of processors, the crack surface evolution in this model was simulated using 512 processors All of simulation s conducted in this study show the robustness of the parallel code in dealing with dynamic fracture problems PAGE 13 13 CHAPTER 1 M OTIVATIONS Fracture is the (local) separation of an object or material into two or more pieces when sufficient work is applied on the a tomic level to break the bonds that hold atoms together (Anderson, 1991) The bond strength is supplied by the interactions between the electrons of atoms. Dynamic fracture problems are those where the inertia forces, rate dependent material behavior, and reflected stress waves can significantly affect the fracture behavior. Since bond rupture at the atomic level is a dynamic process, all fracture is dynamic at the atomic scale. The dynamic fracture problem is the most fundamental in the science of fracture (Cox et al. 2005). The instability of crack propagation is o ne of the most important phenomena in dynamic fracture Experiments have demonstrated unstable crack propagation ( crack branching ) in metallic glass. It is found that the dynamic crack s will osc i llate and branch when they become unstable ( Flores and Dauskardt 1999 ) Although the t heoretical researches have shown that the crack propagation under M ode I loading could become unstable when its speed reaches a critical value (Yoffe 1951, Freund 1974, 1 998, Gao 199 4 ), t he mechanisms for dynamic fracture are still far from well understood. The origin of unstable crack propagation and how a crack will behave when it become s unstable are still open topics currently. Under complex boundary conditions analy tical prediction of the dynamic crack propagation is too difficult to achieve Although the c rack surface under mixed mode loading has been the focus of the community of dynamic fracture for a long time, because of the complexity of the problem, few succes sful works i n this field ha ve been reported Figure 1 1 shows the unstable crack propagation, or crack branching, under PAGE 14 14 the M ode I loading (Flores and Dauskardt 1999) R ecent research (Figure 1 2) shows that the crack propagation will become universal ly un stable under a mix ed mode (I+III) loading (Pons and Kamar 2010). The mechanical origin of this helical crack front instability is still poorly understood. As shown in Figure 1 3, dynamic fracture is a multiscale phenomenon which concurrently involves the b reakdown of material on different length scales. To understand the origin s and mechanism s of dynamic fracture, it is necessary to investigate the failure of material s on different length scales. At the nanoscale or even mesoscale, the atomic structure of m aterial plays an important role in its fracture behavior In single crystal s some planes are natur ally weaker than others. Experiments show that the fracture in single crystal s always happens on these weaker planes. As shown in Figure 1 4, the fracture s urfaces of four different single crystals are not smooth when observed at the micro scale The staircase fracture surface s show the easiest fracture planes in the material. In Figure 1 4 (a) ( Pashley 1960 ) the fracture surfaces along {1 1 1} planes are ca used by the dislocations in a gold thin film, a face c enter cubic (FCC) single crystal In Figure 1 4 (b ) ~ (d) (Schultz et al. 199 4 ) the fracture surfaces are the brittle cleavage planes of three different brittle single crystals, NaCl, CaCO 3 and CaF 2 With the development of large scale computing technique s n umerical methods have bec o me some of the most important and efficient ways to study the mechanisms for fracture During the past few decades a large number of numerical methods have been proposed to simulate fracture. The continuum mechanics based methods, such as: PAGE 15 15 extended finite element method (XFEM), cohesive zone model (CZM), and peridynamics have been widely used to simulate the crack propagation problems without taking into account the atomi stic information of the material As one of the most popular approaches for atomistic simulation molecular dyna mics (MD) simulation is useful in simulating dynamic crack propagation. T here are also many multiscale models which aim to simulate systems with larger size s relative to MD simulation s However, a ll the above methods have their limitations Currently, there are no successful numerical methods for dynamic fracture. I n this work a more prospective method is presented and used for the study of dynam ic fracture problems. PAGE 16 16 Figure 1 1 Crack oscillation and branching under mode I loading (Flores and Dauskardt 1999) Figure 1 2 H elical crack front under mixed mode (I+III) loading (Pons and Kamar 2010) Figure 1 3 Frac t ure at differen t length scales (Markus 2010) PAGE 17 17 (a) gold (b) NaCl (c) CaCo 3 (d) CaF 2 Figure 1 4 Fracture of some single crystals (Schultz et al. 1994) PAGE 18 18 CHAPTER 2 B ACKGROUND State of the art of numerical methods for fracture In this work, a new method for coarse graining (simulating a system using fewer degrees of freedom than they actually have ) the atomistic dynamics of fracture is presented Since there have already been a lot of numerical models for modeling fracture, why do we need another one? To answer this question, i n this section a brief review on the state of the art and the challenges of current numerical methods for dynamic fracture is made Continuum Methods Historically many co ntinuum methods based on continuum mechanics have been successfully applied to study fracture at macro scopic scale. Continuum m echanics is a branch of physical sciences concerned with the deformation and the motion of continuous material media under the in fluence of external effects. External effects that influence bodies appear in the form of forces, displacements and velocities that arise from contact with other bodies, gravitational forces, thermal changes, chemical interactions, electromagnetic effects and other environmental changes. The the ory of continuous media is built upon f our basic laws : (1) the conservation of mass (2) the balance of linear m omentum, (3) the balance of angular momentum, and (4) the conservation of energy. They are universal and are result s from our experience s with the physical world. The constitutive relations are constructed to take the nature of different materials into consideration. These relations s hould not violate the second law of thermo dynamics, which is the famous Cla usius Duh e m inequality in continuum mechanics PAGE 19 19 In the past few decades classic al continuum mechanics has provided most of the theoretical and computational tool for understanding fracture. However, the mathematical framework that has been developed for co ntinuum mechanics is not well suited to simulate such problems as crack propagation since the balance equations of continuum mechanics are partial differential equations, and the presence of the spatial derivatives require continuity of the material. Ther efore, cracks within a material have to be treated as surface boundaries of the material. Many numerical methods have thus been proposed to overcome this barrier in the simulation of fracture. Successful techniques include the extended finite element meth od ( XFEM) (Fries and Belystschko 2006, Huynh and Belytschko, 2009), in which a standa rd displacement based approximation is enriched by incorporating discontinuous fields. T he most significant character of XFEM is that discontinuity can exist within elemen ts. The governing equation for XFEM is just the weak form of the traditional balance equations. To allow the representation of crack discontinuities and voids independent of the mesh, the traditional shape function is modified as follow s : ( 2 1) where are the traditional shape functions, is the displacement at node i, is step function which introduces the discontinuity, are the functions which introduce singularity, and are variables caused by the introduc tion of discontinuity and singularity. This method is static and cannot model dynamic fracture. The crack propagation a nd crack path depend on the empirical criteria which should be given before the computation. Another well recognized method in fracture modeling and simulation is the PAGE 20 20 cohesive zone model (CZM) (Needleman 1987, Camacho and Ortiz 1996, Remmers et al 2003, Pa rk et al. 2009). The advantage of CZM is its ability in model ing arbitrary dynamic crack propagation. By inserting interface elements between continuum finite elements along the potential crack path, the cohesive crack propagations can be modeled. The gove rning equation for this method is as following: ( 2 2) where is the displacement, is the first Piola Kirchhoff stress which is a n energy conjugate of deformation gradient is the surface traction on external boundary is the internal traction on internal boundary and is the separation of the internal boundary. As shown in Fig ure 2 1 two kinds of element, a bulk element and an inter f acial element (cohesive element), are used this model. D iscontinuities are introduced by the failure of the cohesive element It is worth while to mention that the cohesive element is infinite ly thin in its undeformed stage So it is very difficult to define deformation (or strain) for this element. Actually, the constitutive relation for the cohesive element is described as a traction separ ation relation. The constitutive relations for the interfacial element and the bulk element can be expressed in the following forms: ( 2 3) where and are potential energy dens ity stored in bulk element and interfacial element, respectively. It is obvious that the constitutive descriptions of these two element s are different. This inconsistency causes the difference of the mechanical behaviors of these two kinds of element. PAGE 21 21 Some laterdeveloped continuum theories have been trying to extend continuum mechanics to more fields of application, such as microscale. Micromorphic theory (Eringen and Suhubi 1964) is a continuum field theory which considers a material body as a continuous c ollection of deformable particles. It is an attempt to associate the microstructure of a material to its macroscale behaviors under the framework of continuum mechanics. Another interesting attempt, nonlocal continuum theory (Eringen and Edelen 1972, and E ringen 1972a,b), is concern ed with the phy sics of material bodies whose behavior at a point in a body material is influenced by the s tate of other points in the body. It is worthwhile to mention another attempt peridynamics theory, to improve the situatio n through re formulating elasticity theory The p eridynamics theory, proposed and developed by Dr. Silling (Silling, 2000), is a nonlocal continuum theory. Compared with traditional nonlocal continuum mechanics, the primary advantage of the Peridynamics app roach for fracture modeling is that it does not require any supplemental relation to dictate crack initiation, propagation velocity and direction, crack branching, and crack arrest. All of these phenomena emerge as the consequences of the equations of moti on and the constitutive model. Like MD simulation s it does not need any special numerical treatment for modeling of arbitrary crack growth. The ad vantage of p eridynamics in simulating discontinuities has been well demonstrated through simulations of the f ragmentation process of brittle materials (Silling and Askari 2005). However peridynamics simulations are always limited by the accurate determination of material constants. The physics of these parameters and how to measure them in experiments are still challenge s currently. PAGE 22 22 Despite their tremendous success continuum methods have their limitations as well One main limitation is the complete neglect of atomistic mechanisms which are sometimes of interest in fracture problems. An other limitation is that t he constitutive models are all empirically derived from macroscale experiment s For these reasons, continuum mechanics is only suitable for macro and micro scales. Atomi stic Simulation Methods There is little doubt that most of the low energy physics, che mistry, materials science and biology can be explained by the quantum mechanics (QM) of electrons and ions and that, in many cases, the properties and the behavior s of materials derived from the quantum mechanical description of events at the atomic scale (Ghoniem et al. 2003). However, solving the Schrdinger E quations for a large system is very computationally expensive. To determine the properties of an ensemble of atoms larger than that which can be handled by computational QM, molecular dynamics (MD) simulation is one of the most popular ways. With the advent of high performance computers, atomic level molecular dynamics simulation has become a powerful tool to study fracture at the atomic scale (Holian and Ravelo 1996, Zhou et al. 1996, Abraham et al. 1994, Abraham 1996, Abraham and Gao, 2000, Buehler et al. 2003, and Buehler and Gao 2006) The basic idea of MD is using empirical, effective interatomic potential energy functions to eliminate all electronic degrees of freedom T he se potentials are repre sented by function s that depend on the atomic configuration (i.e. relative displacement) and the local environment (i.e. electrons). In general, the total effective potential of the whole system is PAGE 23 23 ( 2 4) w here and correspond to two body interaction s three body interaction s four body interaction s, and N body interaction s respectively. This expression converges relativel y quickly. In practice, only keeping the first 3 terms is sufficient for the desired accura cy The formulations for the effective potential serve as the constitutive equation in molecular dynamics simulation. The dynamic evolution of the system is governed b y classical Newtonian mechanics. F or each atom the equation of motion is given by : ( 2 5) which is derived from the classical Hamiltonian of the whole system: ( 2 6) In ad dition to the unique advantage that MD simulation can provide with regard to the atomic scale mechanism of fracture, it possesses no barrier in modeling and simulation of arbitrary crack growth, branching, deflection, and arrest. Different from continuum m echanics, the governing equation for MD simulation treats discrete atoms in the ensemble, so the variables do not need to be continuous in space (they are continuous in time). In this sense, MD simulation is very suitable for simulating the problems with s pace discontinuities, such as crack propagation. The main limitation of current MD simulation is the sp a cial and temporal size achievable Current state of the art supercomputer MD simulations can handle several billions of atoms, amounting to a volume of less than a cubic micron. There is still a gap between experiments and MD simulations. PAGE 24 24 Multiscale Simulations Methods Motivated by the limitation of MD simulation s tremendous efforts have been made in developing some methods that are both efficient and ac curate. T hese approaches are called multiscale methods. Usually, multiscale methods are classified into sequential methods and concurrent methods. Sequential methods (also called hierarchical methods) attempt to piece together a hierarchy of computational approaches for different length scales. In this methodology, the fine scale model (atomistic model) and coarse scale (coarse grained model) model are weakly coupled and run sequentially. The parameters (or information) obtaine d from the fine scale are pass ed to as an input of the coarse scale. Concurrent methods are those in which the fine scale model is embedded and strongly coupled with the coarse scale model ; t he fine scale model and coarse scale model are always run concurrently. In concurrent simulatio n s the system is often partitioned into domains characterized by different scales. Sequential multiscale simulations : For the applications in different fields, a large number of sequential multiscale models have been developed. Here we only concentrate o n those approaches for fracture problems One of the most significant approaches of sequential multiscale methods aim s to simulate the fracture dynamics, which couple the atomic scale and macroscale, is proposed by Yamakov and colleagues (Yamakov et al. 20 06). In their work crack propagation at the grain boundary of two aluminum crystals was simulated by a cohesive zone model and a traction separation law is extract ed from MD simulation results It is worthwhile to note that this methodolog y is more suitab le for determining the elastic properties of material s However, since the descriptions of fracture at the atomic scale and macroscale are quite different, the transferability of the atomistically PAGE 25 25 fitted decohesion law is very poor Recent research shows t hat cohesive laws of interfaces are not sufficient for modeling fracture problems using the cohesive zone model. (Coffman et al. 2008 ) Concurrent multiscale simulations : Fracture dynamics is one of the most c hallenging problems in material science and soli d mechanics because fracture phenomena are governed by processes occurring over a wide range of length scales. In particular, the physics on different length scales interacts dynamically. Therefore, the concurrent multiscale methods with strongly coupled f ine and coarse scales are adequate for the study of fracture dynamics. This category of methodologies in which certain key regions are treated with fine scale models while most of the domain is treated with a coarse scale model is called partition domain m ethod s In the past decades, a number of concurrent multiscale methods have been proposed. The Macroscopic Atomistic Ab Initio Dynamics (MADD) method proposed by Abraham, Broughton and their colleagues (Abraham et al. 1998, Broughton et al. 1999) aims at linking the length scales ranging from the atomic scale through the microscale. In th eir approach, the crack tip region is treated with a quantum mechanical tight bi nding approximation method; the n molecular dynamics simulation is used for the region arou nd the quantum region ; the region far away from the crack tip is treated via the finite element method in the context of continuum elasticity. H andshaking region s the regions between the atomic and the continuum regions, are used to link each two neighbor ing regions. In handshaking region s the governing equations in both of the two regions contribute to the displacement field. As shown in the gover ning equation for this approach PAGE 26 26 ( 2 7) t he region FE/MD is the handshaking region for FE and MD, the region MD/TB is the handshaking region for MD and TB. Actually, MD and TB regions are fully overlapped which means that the TB region is a subset of the MD region A primary concern in this method is the seamless coupling of force and displacement at the interface between the two regions. The inherent mismatch of the description between different length scales (such as atomistic and continuum scales) can lead to various numerical and theoretical difficulties. For example, the descriptio n of continuum region s and atomic region s are different ; this difference prevents the propagation of high frequency waves from atomic to continuum region s Instead, high frequency waves are reflected and can melt the atomic region An other kind of concurre nt model is based on the Cauchy Born rule or approximation which states that in a crystalline solid subject ed to a small strain the positions of the at oms within the crystal lattice follow the overall strain of the medium. Quasicontinuum (QC) is a famous example in this kind of model QC wa s firstly proposed by Tadmor, Ortiz and Phillips (Tadmor, Ortiz and Phillips 1996) and has been used to model the st atic fracture of nickel (Miller et al. 1998) In the ir QC models an atomic region (fine region) is concurrently linked to a coarse region. E nergy minimization is used for the whole specimen to determine the displacements of the atoms in the fine region as well as the element nodes in the coarse region. In the coarse which represents a group of atoms. When the elem ent size is reduced to an atomic size, there is no difference between PAGE 27 27 nodes and atoms. Motiv ated by this approach, Xiao and Belytsko (Xiao and Belytsko 2004) studied the static fracture of a graphene sheet through a proposed bridging domain method (BDM) which links the molecular mechanics region and the QC region by using Lagranging multipliers. Because they lacking the balance equations for energy, the a fore mentioned approaches as a series of energy minimization methods, can only simulat e static, zero temperature problem s There are two common pitfalls in all kinds of partition domain methods. ( 1) Only static and zero temperature problems can be investigated ( 2 ) The critical phenomenon should be confined in the fine scale region and other parts of the specimen only sever as boundary conditions. So these methods are not suitable for the investiga tion of dynamic fracture Summary In this section, the s tate of the art of numerical methods for fracture is reviewed. From the review, it is seen that the modeling of dynamic fracture is still a challenging topic. A new model which can take into account not only the large scale fracture dynamics but also the atomistic mechanism is required. Most of current models only meet one of the two requirements. The goal of multiscale methods is to overcome this difficulty. However, both sequential and concurrent mu ltiscale methods still have their own difficulties currently. One of the most significant barriers to achiev ing this goal is the inconsistency of the descriptions between atomic and continuous region s The r ecently developed atomistic field theory ( AFT) (C hen and Lee 2005, Chen 2006, 2009) provides a way to achieve this goal In this work a new coarse grained model based on AFT is introduced and used to deal with the dynamic fracture problem Some details of this method will be presented later in Chapter 3 We will also show that PAGE 28 28 the atomistic dynamic s of fracture can be treated accurately with a coarse grained model Introduction to p aralle l computing What i s Parallel Computing? Parallel computing ( from Wikipedia) is a form of computation in which many calculations are carried out simultaneously, operating on the principle that large problems can often be divided into smaller ones, which are then solved concurrently ("in parallel"). If a computational problem has exploitable concurrency, then this problem can be decomposed into several parallel sub problems. These sub pro blems can be executed separately. By using multiple processors to compute these sub problems concurrently, we can obtain the solution faster than only using one processor for the computation of the whole problem Also, if each processor has its own memory, partitioning the data between the processors may allow larger problems to be handled faster than could be handled on a single processor. Why Do We Use Parallel Computing Here? Scientific computations are cost ly in two domains: the computational time and t he memory. These barriers can be overcome by using parallel computing. With multiple processors working together, we can always compute larger systems in shorter time. Two of the successful examples are the applications of parallel computing in molecular d ynamics and finite element method Most of the popular MD computation software such as AMBER, LAMMPS, DL_POLY etc., have been well parallelized. For those well known finite element method packages such as ANSYS, ABAQUS, MSC/NASTRAN et c ., parallel perfo rmance is also an important standard for the evaluation of these products. The better the parallel performance is, the larger the system that can be PAGE 29 29 computed with a fixed amount of processors In the current work, a field theory is implemented by finite el ement method To simulate large systems, a parallelization for the finite element method is conducted. Then the parallel code is used to investigate dynamic crack propagation in crystals. Since d ynamic crack propagation is a multiscale phenomenon which inv olves local material behavior the mesh for the computed model should be fine enough to capture the local information. Since a large amount of memory and computational time is required for the computation of the model, parallel implementation is necessary for the current work. Another reason for using parallel computing here is the high concurrency of the problem. For the current method over 9 5 % of the computations time is taken by the calculation of the internal forces in the finite elements. Using the ex plicit time integrate strategy, the calculation of the internal force on each element node can be solved independently. This means that the step for internal force computation can be parallelized with high efficiency. So the current method is very suitable for parallel computing. High Performance Computing Resource s Currently, there are thousands of high performance computational resou rces all over the world. supercomputer changes about every six months. In June of 2010 the number one supercomputer was JAGUAR located at Oak Ridge National Laboratory in U S. Five months later, TIANHE A1 in China overtook it. By June of 2011 the top ranking supercomputer was changed to K supercomputer in Japan The fast development of s upercomputers indicates that high p erformance computing is being more and more important in research and industry. PAGE 30 30 For our current research, the available high performance computing resource s include the high performance computer (HPC) center in Universi ty of Florida (UF) and TeraGrid. The HPC of UF runs several clusters which is open to the faculties and students of UF. This Center runs several clusters with over 3,000 cores in dual core, quad core, and hex core servers. Most of the servers are part of o ne of two distinct InfiniBand fabrics. The clusters share over 120 terabytes of distributed storage via the Lustre parallel file system. The performance of the HPC in UF is about 15 teraflops. One teraflop means the ability to do 1 x 10 12 floating point oper ations per second. Each dual core, quad core, or hex core server has its own memory from 4GB to 64GB. TeraGrid is an open scientific discovery infrastructure combining leadership class resources at eleven partner sites to create an integrated, persistent c omputational resource. Using high performance network connections, TeraGrid integrates high performance computers, data resources and tools, and high end experimental facilities around the country. Currently, TeraGrid resources include more than 2 petaflop s of computing capability and more than 50 petabytes of online and archival data storage, with rapid access and retrieval over high performance networks. One petaflops means the ability to do 1 x 10 15 floating point operations per second. Researchers can als o access more than 100 discipline specific databases. With this combination of resources, TeraGrid is the world's largest, most comprehensive distributed cyber infrastructure for open scientific research. Message Passing Interface (M PI) During parallel co mputation, it is always necessary for processors to communicate with each other for sharing their computation results dynamically. MPI, a parallel PAGE 31 31 programming environment, provides a convenient way to achieve this goal. MPI created in the early 1990s, is the standard programming environment for distributed memory parallel computers. The standard defines the syntax and semantics of a core of library routines used for writing portable message passing codes in Fortran or C programming language s Currently, there are several implementations of MPI such as OpenMPI, MPICH, HPMPI, IntelMPI, and Microsoft MPI. Among these implementations, OpenMPI and MPICH are the most popular. Each implementation includes some compilers which are used to compile the source code in different language s For example, in OpenMPI, mpicc is the compiler for the c ode written with C; mpif90 is the compiler for the code written in Fortran 90; mpif90c++ is the compiler for the code wri tten in C++. A n implementation is a package installed in a system and a compiler is used for compiling a source code written by the use r. The executable generated by the compiler can be run There are s ome important concepts introduced in the first version of MPI, such as environment management, point to point communication, collective communica tions and communicator. For the second version, MPI 2, some extensions, including parallel I/O, one side operation, and d ynamic process management, are provided These concepts are closely related to the application of the MPI standard. A portable code wri tten with MPI must be compatible for different systems. To achieve this portability, a subroutine called MPI_INI T ( ), which is used for the setup of the parallel environment need s to be called before any other MPI routines. At the end of the parallel code another subroutine, MPI_FINALIZE( ), is used to clean up all MPI PAGE 32 32 states. Point to point communication, the transmittal of data from one process to another, is the basic communication mechanism of MPI. In this type of communication, a SEND operation is al ways used combined with a RECV operation. There are two kinds of point to point communication modes: blocking and unblocking. For blocking mode, the SEND and RECV operation s are always finished at the same time. For unblocking mode, one of these two operat ions could be finished first. Then a subroutine, MPI_WAIT( ), is used before the next communication to avoid conflict in data transmittal Collective communications provide a convenient way to transmit data among a group of processes. The functions of coll ective communications include: synchronization, broadcast ing gather ing scatter ing, and reduction. From our own experience s it is found that the collective communications always take a long time when the number of processes is large. A communicator is an object that specifies a communication domain for process groups. It includes the management of groups, the communications within or between groups. At program startup, the runtime system creates a common communicator called MPI_COMM_WORLD. MPI 2 provides subroutines for transferring data to or from files on an external storage in a high performance parallel way. The files can be randomly or sequentially accessed using the parallel I/O. For the sequential mode, the file pointer for read and write should be Regular send/receive communication requires matching operations by both the sender and the receiver. Remote Memory Access (RMA) an extension in MPI 2, allows PAGE 33 33 one process to spe cify all communicat ion parameters, both for the sending side and for the receiving side. This communication mechanism provides more flexibility for the parallel coding. The MPI 2 process model allows for the creation and cooperative termination of processes after a n MPI appl ication has started. A process can start new processes and establish communication for them by calling the subroutines MPI_COMM_SPAWN ( ) and MPI_COMM_SPAWN_MULTIPLE ( ) The new processes will be cooperatively terminated if the subroutine MPI_INIT( ) is cal led. Running a Parallel Code on High Performance Computers R un ning a parallel code may be different depending on the systems. However, all the systems have a similar structure which is composed of submit nodes, computing nodes, and visualization nodes. Su bmit nodes are for submitting jobs, compiling codes, and file transferring. Computing nodes, being the largest part of the system, are used for running jobs. Visualization nodes are often used to visualize the results wh ich could contain an extremely larg e amount of data After an executable is generated, a script file is used for submitting the computation job. In general t he script file should specify the characters of the job, such as: the name and the path of the executable, the path for the I/O, the number of cores to be used, and the memory required. Then the job is submitted using the command scrip file name Summary With the fast development of high performance computing technology, parallel computing is being widely applied in the field of scientific computation This technology allows us to investigate the real mechanisms more accurately and reasonably with the PAGE 34 34 aid of com puter MPI is a parallel programming environment that provides a convenient way for the implementation of parallel co de. Our current high performance computing resourc e s include the HPC in UF and the TeraGrid. In this work, a high performance parallel code is developed using the program ming language FORTRON in the MPI environment. Figure 2 1 Cohesive element and bul k element PAGE 35 35 CHAPTER 3 A COARSE GRAINED MET HODOLOGY BASED ON TH E ATOMISTIC FIELD THEORY (AFT) Introduction In this chapter a newly developed predictive coarse grained method based on the Atomistic Field Theory ( AFT) (Chen and Lee 2005; Chen 2009) is introd uced AFT is a continuous field theory deri ved solely from atomistic model. This theory provides a new framework for coarse grained and multiscale simulations. Within this framework, both fine and coarse regions are treate d by the same set of governing e quations which will be introduced in some details later B ecause of the consistency of the description, the transit ion from fine to coarse regions becomes natural. The finite element implementation of AFT leads to a new CG method The key feature that ma kes this method different from the conventional finite element method is the construction of the finite element For the new CG method based on AFT, the finite element mesh design is closely related to the lattice structure of the material. Two different f inite elements corresponding to 2D triangular lattice and 3D face center cubic (FCC) crystal are introduced. The internal force field is directly derived from the interatomic potential. Both node integration and Gauss integration strategies for the calcula tion of nodal forc e are presented in this chapter A brief review of coarse grained (CG) models In general, CG models indicate those models with fewer degrees of freedom than they actually have (Espanol 2004). It is well accepted that MD simulation s can pr edict the properties and behaviors of various systems at atomic level from first principals. However, the limitation of MD is also obvious. It is limited by the special and temporal size of the simulation system. So CG models are necessary for the systems larger than PAGE 36 36 that can be handled by MD simulation One of the most intuitive way s to construct a CG model is to treat a group of atoms as a rigid particle. Then an effective force field for the CG particles is used instead of the interatomic potentials. By grouping atoms together, the total number of degrees of freedom can be significantly reduced. There are two key concerns that affect the accuracy of a CG model: (1) the selection of the center of CG par ticle, and (2) the construction of the effective force field. There are many ways, or criteria, proposed to deal with these two concerns For example, minimizing the error of distributions and thermodynamic properties (Soper 1996, 1997, Lyubartsev and Laaksonen 1995, 1999); or using the force matching (FM) ap proaches (Izvekov et al. 2004, Izvekov and Voth 2005, Hone et al. 2005). These models have been widely used in the field of macro molecular soft materials where some of the interatomic bonds are much stronger than others. However, since the deformation of each CG particle is constrained this class of models is not suitable for crystalline material where all the bonds are at the same level. B ase d on the assumption of Cauchy Born rule, there is another class of CG methods which is suitable for crystalline ma terial. In these models, only a small portion of atoms in the system are selected to represent the whole system. The selected atoms are called repatoms (representative atoms) This repatoms determine the displacement field in the system. The motions of oth er atoms are approximated by some collective motions. As is mentioned in C hapter 2 qusicontinuum (QC) method falls in the scope of this category. All of the CG models reviewed a bove are particle based models whose theoretical framework is discrete particle dynamics. There is another class of CG models based on PAGE 37 37 field representations. Classical continuum mechanics as one of the earliest coarse grained field theories, view a material as a homogeneous and continuous medium. The basic or the smallest st ructural unit is idealized as a point mass without structure. All the quantities in the system can be described by field functions. These density functions are related to each other by governing equations. The expression s for these functions can be found b y solving the governing equations analytically or numerically. In continuum mechanics, t he atomistic information, such as the atomic structure and the internal deformations of materials are ignored. So the description of material s is based on the macrosco pic, homogeneous motion of the material particles. For this reason, the material properties involved in continuum mechanics cannot be link ed to the motion of atoms at nanoscale. The most common way to obtain these material properties is from macroscopic ex periment s In other words, this model is not predictive. The fundamental goal of statistical mechanics is to link the microscopic many body dynamics to the phenomenological averaged continuum description. Equilibrium classical statistic mechanics provides a successful link between the atomic motion and the static, equilibrium macroscopic variables (Gibbs 1902). The nonequ ilibrium statistical mechanics extends this link to dynamic situation s The most significant works in the field of nonequilibrium statisti cal mechanics are a series of articles of Kirkwood and his co operators (Kirkwood 1946, 1947, Kirkwood et al. 1949, Irving and Kirkwood 1950, Irving and Zwanzig 1951, Zwanzig et al. 1953, Zwanzig et al. 1954, Ross and Kirkwood 1954). In their works, the equ ation of continuity, the equation of motion, and the equation of energy transport are expressed in terms of molecular variables. However, classical statistical mechanics has only worked for macroscopic properties PAGE 38 38 from the simplest lattice systems work in 1982 (Hardy 1982) i s an e xtension of theory. With the new formulas, some difficulties arising in the old theory can be avoided They make the old theory easier to use in the practice. An introduction to AFT Decompositi on of the Atomic Motion Atomistic Field Theory is a generalized continuum theory. Different from traditional continuum mechanics, in AFT the atomic information is associated with each material point. Figure 3 1 shows the relationship between traditional co ntinuum mechanics, micromorphic theory and AFT. In classical continuum mechanics, each material point is treated as an infinitesimal point, the degree of freedom for this point is 3. In the later developed micromorphic theory, the material point is not lon ger an infinitesimal rigid point but a deformable particle. This particle has nine independent degrees of freedom describing both stretches and rotations, in addition to the three classical translation degrees of freedom. The recent developed AFT views a c rystalline material as a continuous collection of lattice points, while embedded within each point is a group of discrete atoms. Thus the degree of freedom for each lattice point depends on the number of atoms associated to this point. It is worthwhile to mention that the governing equations for AFT could be based on time averaged quantities or instantaneous quantities. The governing equations based on time averaged quantities are derived by Chen and her cooperators in 2005 (Chen and Lee 2005) and the gover ning equations based on instantaneous quantities are derived by Chen in 2009 (Chen 2009). It is also interesting to find that both the time averaged quantities and the instantaneous quantities are following the same set of governing equations. In this diss ertation, the instantaneous version is chosen for the PAGE 39 39 i ntroduction of AFT. One of the key issues of AFT is the decomposition strategy which allows the concurrent treatment of the crystalline material at different length scales. As shown in Figure 3 2, the atomic position is decomposed into lattice position and the relative internal position : (3 1) where is the position of the atom in the unit cell, is the position of the unit cell, and is the relative position between the atom and the centre of the unit cell. The correspondence of this decomposition in the physical space is: (3 2) where z, x, y are the counterparts of and respectively. The Local Density and Its Time Derivative The link between the local density function in physical space and any measurable dynamics function in phase space is established as (3 3) where N is the number of lattice cells and n is the number of atoms within a lattice cell; or is the local density function corresponds to ; is the localization function which can be a Dirac function or a distribution of weighting function and has a unit of inverse volume. Note that the function is continuous for both x and y. The cell averaged continuously distributed local density function is defined as: PAGE 40 40 (3 4) In the above equation, V( x ) is the volume of the lattice cell that is located at spatial place x The condition is derived from the normalization conditions : (3 5) Let represent the internal relative position of th atom embedded within t he lattice cell. It is obvious that the volum e of a unit cell is the sum of the volume of each atom within the unit cell. So Equation (3 4) can be further written into the following form (3 6) From the above equation, a new atomic scale local density is defined as (3 7) This new local density function is continuous in x while discontinuous in The time evolution of this function can be obtained in the spattial description directly from its PAGE 41 41 definition as (3 8) where and From the definition of localization function, one can readily prove that (3 9) Then it follows that (3 10) Bal ance Equations Following from E quation (3 7), the continuous local mass density the linear momentum density and the total energy density in a multiatom system can be defined as PAGE 42 42 (3 11) (3 12) (3 13) where is the site potential energy and is the atomic level local velocity field. With the abo ve three expressions, the conservation of mass, the balance of linear momentum and the conversation of energy equations can be derived following the time evolution laws given in E quation (3 10). It is also proved that the balance of moment of momentum equa tion is automatically valid in AFT (Chen and Lee 2005). Conservation of mass: (3 14) Balance of linear momentum: (3 15) Conservation of energy: (3 16) where and are the homogeneous and inhomogeneous parts of moment flux; and are the homogeneous and inhomogeneous parts of heat flux, respectively. They are defined as (3 17) PAGE 43 43 (3 18) (3 19) (3 20) where and are the difference bet ween phase space velocities and physical space velocities; for pair potential is the internal force acted on atom ; is the external force; and (3 2 1) From E quation (3 17) ~ (3 20), it is seen that and are respectively the momentum flux and heat flux due to the motion and deformation of the lattice. They are consistent with the stress (Cauchy stress) and heat flux concepts in the mechanics of homogenized continua. and are respectively the momentum flux and heat flux due to the reorganization of atoms within the lattice cells. They describe the contr ibution to stress and heat flux due to atomic scale inhomogeneity. An instantaneous expression for the temperature in classical many body dynamics is expressed in terms of the mean squared velocity relative to the local velocity field (Hoover, 1986,1991) (3 22) PAGE 44 44 where is the volume of the unit cell located at x is the total mass of the atoms in the unit cell, is Boltzmann constant With this tempera ture definition and also the internal force density the balance equation of linear momentum ( E quation (3 15)) can be rewritten as (3 23) or (3 24) where is the displacemen t of the th atom embedded at the material point at x ; is the fraction of contribution of the th atom to the temperature. Generally, the internal force density (3 25) is a nonlocal and nonlinear function of atomic displacements, and can be obtained through fitting to experimental measurements or first principle calculations. Equation (3 25) serves as the constitutive equation which links the displac ement field to the internal force. It is seen from E quation (3 25) that the internal force may depend on the high order gradient of displacement field which corresponds to the strain and strain gradient in classical continuum mechanics. Recall that the E q uation (3 14) ~ (3 16) have the similar form as the conservation of mass, balance of linear momentum and conservation of energy equations in classical continuum mechanics, respectively. For this reason, we can treat AFT as a generalized continuum theory. D ifference between Particles Dynamics and Continuous Field Theories Unlike many body dynamics, a continuous field theory is only concerned with local PAGE 45 45 densities, such as mass density, stress tensor, heat flux vector and energy density. The evolutions of l oca l densities in space and time are governed by balance equations. Supplemented by constitutive relations, the balance equations completely determine the dynamic behavior of a material system under a given external field and properly imposed b oundary and ini tial conditions. With the local densities defined in the continuous field theory, balance equations, always serve as the governing equations, can be derived and used to describe the motion of the system. While for particle dynamics, classical mechanics ser ves the governing equation for many body discrete systems. AFT is a continuous field theory which is based on local densities and balance equations. Different from atomistic simulation methods which are based on discrete many body dynamics, temperature is defined as an independent variable in AFT. This variable can be directly solved from the conservation of energy equation. It is interested to find that the balance equations in AFT will be equivalent to many body dynamics for MD simulation when focusing on the motion of every atom in the system. The finite element implementation of AFT Weak Form of the Governing Equation From the previous section, it is seen that AFT is a nonlinear and nonlocal field theory which is able to describe the motion of materials at any length scale. The governing equations of AFT have been numerically solved by meshfree method to investigate the dynamic fracture of a model composite material (Deng et al. 2009) In this dissertation, a finite element implementation of AFT is introd uced. It is seen that, for materials under a given temperature field, Equation ( 3 24 ) supplemented by Equation ( 3 25 ) completely govern the evolution of the atomic scale displacements in space and PAGE 46 46 time. Substituting E quation (3 25) into E quation (3 24), we have the following governing equation (3 26) here we omit the temperature term since we are considering a homogeneous or a given temperature field. Under these special situations, the temperature term which has been omitted h ere is either zero or a constant. Therefore only E quation (3 26) is needed as the governing equation. The conservation of mass equation and the conservation of energy equation (Equation (3 14), (3 16)) are automatically satisfied. There are 3n equations in (3 26) need to be solve d for determining exactly the same number of unknowns ( =1, 2, n). In this section, the finite element method is used to numerically solve E quation (3 26). Notice that E quation (3 26) has the similar form as the governing equation for classical continuum theory. Different from the classical continuum theory which explicitly involves the spatial derivatives of displacement field in evaluating the internal force field, E quation (3 26) employs the integration over the domain covers the whole system, to determine the internal force fie ld. Benefit from the avoiding the spatial derivatives, we can easily overcome the difficulties in solving the high gradient or local discontinuity problems with AFT. F irst ly, we approximate th e displacement field within each element by (3 27) whe re stands for the approximation of atomic displacement field within the element, is the shape function, and is the i th component of the displacement of th e atom within the node. Here, as shown in Figure 3 3, each element PAGE 47 47 node is associated with a unit cell which contains n atoms. Then t he weak form of E quation (3 26) can then be written as (3 28) where both and represent the domain of the whole system. And on the left hand side of this equation, the first term is the inertial term; the second term is the internal force term; the third term is the external force term. Numerical Evaluation of Each Term in Weak Form Equation To calculate the inertial term, the lump mass assumption is used. T he lumped mass matrix can be defined as the following form: (3 2 9) The internal force term can be obtained with numerical integration strategies such as Gauss integration or nodal integration. For nodal integration, the discretized internal force term can be written in the following form: ( 3 30) where i s the summation over all the nodes; is the weight associated with each node. It is seen from the above equation that the internal force term is a function of the displacements of all the nod es. PAGE 48 48 It is easy to see that the integration of the internal force over each element is contributed by two different parts: internal force on the boundary and in the bulk region within each element. The reason for this difference is that the internal force in the bulk region of an element only depends on the nodal displacement of this element, while the internal force on the boundary of an element not only depends on the nodal displacements of this element but also depends on the nodal displacements of its n eighboring elements as well. In other word, the bulk region part only locally depends on the displacement field within the element, while the boundary part also depends on the displacement field of the neighboring elements. For homogeneous deformation the displacement field varies smoothly with space, so there are few differences between interior and boundary regions However in some critical situation, such as dislocation or crack nucleation and propagation, the displacement field may have high gradient i n some critical regions. At this time, it is necessary to distinguish the interior and boundary regions (3 31) T he above E quation (3 31) is for the numerically evaluation of internal force using Gauss integration strategy T he integration of internal force over the whole domain is divided into three parts: the nodal part, the boundary part and the bulk region part. The node part includes all the internal forces on all the element nodes in the system. PAGE 49 49 means the s ummation over all element nodes. The boundary part includes the internal forces on the boundaries of all elements. means the summation over all the gauss points on the surface. The bulk region part includes the interna l forces inside each element, the internal force field in this region should only depend on the nodal displacement of the element where the region located. means the summation ove r all the gauss points in this region The notatio ns and represent the weights for the surface gauss points and bulk gauss points respectively. It is worthwhile to mention that the weight for the nodal part summation is 1 instead of in nodal integration version. Compared with the Gauss integration procedure, nodal integration is significant more efficient, easy to implement, but less accurate and usua lly leads to instability In order to stabilize the nodal integration strate gy, the nodal integration form can be rewritten, into the following form (3 32) In the above Equation (3 32), since the integration domains and both indicate the region that cover the whole system, so there is no difference to exchange the indexes in tensor PAGE 50 50 analysis. So Equation (3 32) can be further rewri tten as (3 33) In summary, Equation (3 33) is employed to calculate the internal force on the th atoms of the th node in nodal integration version. While in the Gauss integration version, Equation (3 31) is used for the com putation of internal force. Figure 3 4 shows the differences of these two numerical integration strategies. In the pictures, the solid dark points indicate the nodal points where both the internal forces and the displacement nee d to be determined and the s tar points indicate the Gauss points where only the internal forces need to be determined. Since the external force field is always explicitly given, so the external force term can be analytical ly derived out. It is also worthwhile to mention that the tem perature term is omitted in this work because of the temperature field is always assumed to be given in our current consideration. If the temperature field is not explicitly given, the thermal mechanical coupling problem needs to be solved. At this time, t he equat ion of conservation of energy (Equation (3 16)) will be involved to solve the additional field PAGE 51 51 variables, temperature. Finally, the discretized gover ning equation (3 34) is derived by combining Equation (3 28) ~ (3 33) The central difference method is used for sol ving the discretized governing E quation ( 3 33 ). Some details of the central difference algorithm in solving this equation will be addressed later in C hapter 4 Mesh Strategies Since there is no space de rivatives involved in the formulation, the continuum finite elements do not have to be connected. This provides us with more freedom in designing the mesh strategy. As shown in Figure 3 5, there are tw o options for the mesh strategy: with and without conne ctivity. Of course the latter is suitable for simulating with coarse mesh, the nucleation and propagation of discontinuities throughout the whole region. But since the coarse mesh is used instead of the finest mesh, some of the high frequency waves might be suppressed when using low order interpolation functions. This unphysical suppressing can be avoided by using physically based shape functions. In this work, our goal is to simulate the artificial crack propagation in material. So the coarse mesh without connectivity strategy is employed here And since the element size is very small, only linear shape functions are used. It is worthwhile to mention that since there is no connectivity between neighboring elements, each node only belongs to one element. Th en the integration over the whole domain can be divided into the integrations over the elements in the system. Thus the summation of all nodes or Gauss points in the system can be simply divided into the summation of the nodes or Gauss points in each eleme nt. In other word, no assembling is needed. PAGE 52 52 Since the single crystalline material is anisotropic, its fracture behaviors are direction depended. For these reason those planes which are easier to be separated should be chosen as the surfaces of element. Di fferent from in the cohesive zone model, these potential paths for fracture evolutions is physically selected. This is also one of the significant advantages of the new CG method over the classical finite element method. In the following, two different fin ite elements based on material structures are introduced. 2D triangular element In the study of 2D brittle fracture problems, 2D triangular lattice is the most frequently used model (Holian and Ravelo 1995, Zhou et al. 1996, Abraham et al. 1994, Buehler et al. 2003, Buehler and Gao 2006). This lattice model is simple and useful. The model represents the (1 1 1) plane of the face center cubic (FCC) crystalline material. Since the triangular unit formed by the neighboring lattice points is equilateral triangl e this lattice model is isotropic. For a triangular lattice shown in the left side picture of Figure 3 6, the density of bonds along the red dash lines are less than the density of bonds along the blue dot dash line. This means that the crack is easier to propagate along the red dash lines than the blue dot dash line. Actually, the bond densities along these three directions are the lowest in all the possible directions. So these red lines outline the three possible directions of weakest planes. The right side picture of Figure 3 6 shows the shape of a finite element. This element is in the shape of equilateral triangle to maintain the isotropy of the material. As shown in Figure 3 6, area coordinates are used as the shape functions for this 2D triangular e lement. The red stars in the triangular element indicate the Gauss points for edge and area integrations. Then following the traditional procedure for treating the PAGE 53 53 triangular element, the shape functions are substituted into Equation (3 34) to get the equa tions for solving the nodal displacements. 3D hybrid element It is widely accepted that FCC crystal s seldom show brittle behaviors. In most of the case, this kind of material s display ductile properties with dislocation dominant its fracture. Zhou et al. ( 1999), through MD simulation, showed that in FCC crystal the crack propagation is close related to the dislocation. They also obs erved that the crack always kink s to the directions along the slip planes. Motivated by this finding, a hybrid element composed of two tetrahedrons and one octahedron is introduced here to study the fracture of FCC crystal. The thick red lines and blue lines in the left picture of Figure 3 7 outline a primitive cell of the FCC crystal. a1, a2, and a3 are three lattice vectors. As shown in Figure 3 7, this primitive cell can be divided into two tetrahedrons and one octahedron. For FCC crystal, every slip planes should be parallel to one of the four surfaces of the t etrahedron shows in Figure 3 7. Obviously it is too cost to embed al l slip planes in a CG model. However, all these slip planes can be sorted by their orientation into four classes. Within each class, the slip planes are parallel to each other. The tetrahedral element in Figure 3 8 shows the four classes of planes. Every s lip plane in this material should be parallel to one of the four faces of the tetrahedral element. As shown in Figure 3 8, a hybrid element for FCC single crystals is composed of two tetrahedral sub elements and two pyramid sub elements. The shape of this hybrid element is the same as a rhombohedral primitive unit cell of the FCC single crystal. By introducing two internal surfaces between the tetrahedral and the pyramid elements, the PAGE 54 54 hybrid finite element, with 16 faces, shows all of the four classes of s lip planes. Note that the 16 faces are divided into four groups and each group, with four parallel faces, corresponds to a class of slip planes. Furthermore the total areas of each group of faces are the same. So the discontinuity can propagate in all the possible directions with equal possibility. Within each hybrid element, the conventional linear shape functions are used to interpolate the displacement field in each sub element. Between neighboring hybrid elements, there is a gap of ( is the length of a lattice constant) to separate different elements as well as the gap separates tetrahedral and pyramid sub elements. Since the displacement field across the gap is not necessary to be continuous, the dis continuities, such as dislocations or cracks, can propagate along the space within or between elements. The shape functions for (3 35) where is the i th component of the volume coordinate at point x Summary In this chapter, a new coarse grained method based o n AFT is presented in details. One of t he most significant feature s of the new developed method is its outstanding ability in simulating dynamic fracture problems without losing the atomic inform ation. As a field theory, AFT unifies the governing equations in both atomic and continuum regions. Under the framework of AFT, the interatomic potential can be directly served as the constitutive relations for solving the balance equations. The coarse gra ined method developed in this work is a finite element implementation of AFT. Instead of stress strain PAGE 55 55 relation, in the new coarse grained method, interatomic potential is used to calculate the internal nodal forces. It is worthwhile to mention that not on ly the interatomic potential but also the atomic structure of the material are embedded in the coarse grained model Each finite element should correspond to a specified material In this chapter, two kinds of finite elements are developed for two differen t materials. A 2D triangular element is designed for 2D triangular lattice. A 3D hybrid element is designed for 3D single crystal with rhombohedral primitive cell. One of the advantages of this treatment is that the faces of the element are always the weak est planes where the failure is most likely to occur on them. This CG method has been successfully used to investigate the dislocation and dynam ic crack branching in material (Xiong et al. 2011; Deng et al. 2010) In C hapter 5 and 6, some details on the ap plication of this CG method for dynamic fracture will be presented. It is worthwhile to mention that, for the present method, only the balance equation of linear momentum is solved since only given temperature cases are under the consideration here When c rack is propagating, the temperature field near the crack tip could become very complex. To take into account the thermo mechanical effect we need to solve the conservation of energy equation along with the balance equation of linear momentum. PAGE 56 56 Figure 3 1 Relationship between AFT and continuum mechanics Figure 3 2 Decomposition of atomic positions Figure 3 3 Schematic scratch for finite elements in the FE implementation of AFT PAGE 57 57 (a) Nodal integration (b) Gauss integration Figure 3 4 Nodal integration and Gauss integration (a) with connectivity (b) without connectivity Figure 3 5 Two options of the mesh strategy Figure 3 6 The triangular el ement which shows the weakest planes PAGE 58 58 Figure 3 7 The p rimitive unit cell and slip planes of FCC crystal Figure 3 8 Sketch of one hybrid element PAGE 59 59 CHAPTER 4 PARALLEL IMPLEMENTAT ION OF CG AFT METHOD Introduction The fast development of high performanc e computers has provided us a great opportunity for massive computation. Using multi processor platforms, not only the computational time but also the memory can be largely reduced These two advantages are also critical to the scale up of the simulation u sing the new CG method. For most supercomputers, the memory allocated to each core is less then 8G B So the total number of elements that can be handled by a series code should not exceed 4000. W ith the parallel code, the total number of elements in a simu lation model can be much larger than this critical value. In this chapter, the structure and the implementation of the parallel code is introduced in details. The performance of this parallel code is tested via stability and scalability analys e s. To test t he stability of the code, the same problem is computed with different number s of processors. For a stable simulation, these results should be identical to each other. The scalability of a parallel code reflects how the parallel code performs with different number s of processors. In this work, the effects of the size of the simulation system on the performance of the parallel code are also investigated. Four design spaces for the parallel code T he re are four design spaces in the design of a parallel code: f inding concurrency, algorithm structure, supporting structures, and implementation mechanisms In this section, some details of these four design spaces and their specification s for the studied code will be presented. In the finding concurrency design spa ce, the tasks and data are decomposed first. PAGE 60 60 Then the decomposed sub tasks are grouped and ordered according to their functions and sequences Finally, the data sharing for the sub tasks is analyzed. This design space, including decomposition and dependenc y analysis is an important step in implementing a parallel program. After decomposing and grouping the key tasks, we have the following five groups of tasks: 1. a group of tasks to read input data 2. a group of tasks to construct the neighbor list 3. a group of ta sks to calculate the nodal forces 4. a group of tasks to update the position and velocities of each node 5. a group of tasks to output the results There are three groups of key data: 1. t he coordinates and velocities for each node 2. t he forces for each node 3. t he neig hbor list for each node D ata group 1 is initiated by task group 1. Then it is used by task group s 2 and 3. Task group 4 update s the data in data group 1 Data group 2 is generated by task group 2 and then is used by task group 4 in the updating Data grou p 3 is generated or updated by task group 2 and then used by task group 3. The data groups used by task group 5 depend on what information need s to be output. The goal of the second design space, algorithm structure, is to decide which pattern or patterns are most appropriate for the problem. For different problem s different algorithm structure pattern s are needed. The criteria for choosing the pattern is to best exploit the concurrency. The program is composed of a large loop of time evolutions. The resul t of each time step is an input for the next step. Within each time step, the nodal forces and positions are updated independently. So we can use different processors to handle these jobs concurrently. For this problem, either task parallelism or PAGE 61 61 a geometr ic decomposition algorithm could be used. The former is easy to implement but need s more time in communication between processors; the latter is more difficult in implementation but could save more communication time especially for very large system s I n this work, to enable our program to handl e an extremely large model, we use the geometric decomposition algorithm. The third design space, supporting structures, is more close to the programming. In this space, the program structure and data structure ar e determined. In this work, the SPMD (single program, multiple data) pattern is used for the tasks within one time step. The distributed array pattern is used to partition the large arrays between different processors. Since the geometric decomposition alg orithm is used, the arrays are partitioned according to the geometry of the nodes. The fourth space, implementation mechanisms, includes the details of the the parallel program implementing The three core subset s of these mechanisms are : process or thread management, synchronization, and communication. As w as mentioned before, MPI is used for the parallel program implementation. MPI is fundamentally based on processes. A n MPI program with an executable named foo is launched through the command: np The processes are destroyed when the programs running on the nodes of the parallel computer exit. T here are a large number of operations in MPI, the most useful of which are in the following three groups: (1) b asic point to point message passing; (2) c ollective operations; (3) a dvanced point to point m e ssage passing. The details of these operations will not be listed here because they can be easily found in any MPI manual such as: MPI The Complete Reference by Marc Snir et al. PAGE 62 62 After the above analy sis, the parallel strategy for this problem is given by the following steps: 1. Read the information from the input file. This should be done by every processor at the beginning of the computation. This is done by subroutine INPUT 2. Divide the whole domain int o P sub domains (P is the number of processors) and find neighbors for each sub domain. Instead of restoring the information for every node in the system, each processor only restores the nodal information, such as positions, velocities and forces, of its own sub domain. Each processor also maintains a local list of the local indices and global indices of the nodes and elements in this sub domain. This is done by subroutine GETMAP 3. For every few time steps, construct the neighbor list with the link cell str ategy. Each processor only constructs the neighbor list of the nodes in its own sub domain. The radius of this neighbor list is This is done by subroutines ELIST 4. For each time step, update the positions of nodes of boundary ele ments by inter processor communications. This is done by subroutine EXCHANGE_U 5. For e ach time step, compute the forces of each node within the sub domain. This is done by subroutine GETFORCE 6. For each time step, update the positions of nodes within the sub domain. This is done by subroutine CENTRAL 7. Output the results every N steps ; N is set by the user of the program. This is done by subroutine OUTPUT The structure of the parallel code The flow chart shown in Figure 4 1 depicts the structure of the curren t parallel code. In this figure, t_GMAP, t_ELIST, t_EXC, t_CFR, and t_CEN represent the CPU time for the key steps GETMAP, ELIST, EXCHANGE_U, CFORCE, and CENTRAL, respectively. EXCHANGE_U, CFORCE and CENTRAL are the three steps which need to be executed ev ery time step. The average number of time steps for a calculation is on the order of several tens of thousands. The time t_SOLVE equals to the sum of t_EXC, t_CFR, and t_CEN. In the following, these key steps will be introduced separately. PAGE 63 63 Domain Decomposi tion GETMAP GETMAP is used to divide the whole system into several domains w it h no overlap s Each domain is allocated to a single processor. The more evenly assigned the load is among the processors, the more efficient the parallel computing. Therefore, i n this subroutine the whole system is divided into domains with almost the same number of finite elements. To divide the model shown in Fig ure 4 2 (a), first a mapping protocol is used to transfer the model in to the cubic shown in Fig ure 4 2 (b). T h en the cubic is evenly divided into many domains as shown in Fig ure 4 2(c) In current CG AFT method, solving the nodal positions of each domain involves not only the nodes within this domain but also the nodes of its neighboring domains. Figure 4 3 shows the interfaces between domains T he interactions across these interfaces involve the neighboring domains. Thus the other objective of this subroutine is to record the neighbors of each domain which will be used in the future for data communication between nei ghboring domains. Construct Neighborlist, ELIST As is presented in C hapter 3 AFT is a nonlocal theory and the internal force on each atom de pends on its neighbors through E quation ( 3 25 ). In the proposed coarse grained method, only the forces on Gauss po ints are needed for numerically evaluating the integration. I t is therefore necessary to construct a list of neighbors, called neighborlist for each Gauss point before the calculation of its internal force. ELIST is a subroutine designed for the constru ction of neighborlist. Following the numerical integration scheme introduced in S ection 3. 2, there are Gauss points associated with node, edge, surface, and body. ELIST will construct the neighborlist for all the Gauss points. ELIST is divided into two mai n steps T he first step is to construct the element PAGE 64 64 neighborlist for each element T he second step is to find the neighborlist of the Gauss points of each element only for its neighboring elements. Here the assumption is made that the finite element is lar ge enough so that all the neighbors of a Gauss point are locate d in its own element or the neighboring elements. Considering a system with m finite elements, if there are n unit cells and g Gauss points in each element, then the computation complexity of E LIST is of the order Here the factor 27 indicates th at each element has 27 neighbor ing elements in a 3D system For a non uniform mesh, this number might vary slightly For a parallel scheme, each processor only deals with the Gauss points within its own domain. If there are n o processors used, then the computation complexity for each processor is of the order When m is very large, the term is much smaller than the number 1 Then t his ste p involves order calculation s However, it is not necessary to call this subroutine every time step. In most case s it is called every few hundred or thousand steps. Communication Between Processors, EXCHANGE U EXCHANGE_U and GETMAP are two of the most important steps in parallel implementation. GETMAP provides the parameters for communication and EXCHANGE_U performs the communication with the parameters provided by GETMAP. As shown in Fig ure 4 4, EXCHANGE_U communicates the updated nodal displacement s between neighboring domains. Fig ure 4 4 shows an example of 3D communication. This communication is divided into three steps and in each step the PAGE 65 65 communication is performed in one dimension (x direction, y direction or z direction) A fter these three step s the information of each domain is shared by all of its neighboring domains. Calculat ion of Nodal Forces, CFORCE The step of calculati ng the nodal force, CFORCE, is the implementation of numerical int egration Using E quation ( 3 31 ) o r (3 33) the Gauss integration or n odal integration strategies are used to calculate the nodal force numerically In this step, th e nodal force of a finite element is calculated based on the internal force s on all Ga uss points belong ing to this finite ele ment. Each processor only needs to handle the elements assigned to it T his is the reason why the computational time and memory can be largely saved by using the proposed parallel algorithm The efficiency of this step determines the overall efficiency of the parallel code. This will be discussed in detail later in the following sections. Time Integrations, CENTRAL The introduced new CG method involves both spacial and t ime integrations The spacial integration has been discussed above T he time integration is numerically achieved using CENTRAL. In this work, the explicit numerical integration scheme, the central diffe rence method, is applied for th is time integration. Like in CFORCE, each processor only needs to handle the nodes within its own domain. Fig ur e 4 5 shows the main steps for the central difference method used here. Performance of the parallel code Stability Stability is one of the most important features of a computational code. For parallel code, the computational results are always different w ith different number s of PAGE 66 66 processors. The computation results from an unstable code are useless. Unstable computation always results in nonphysical or infinite ly large values s o the stability test is the first step in test ing the performance of our paralle l code. As the first example, the stress wave propagation in a face center cubic (FCC) single crystal bar is simulated. As shown in Figure 4 6, 6400 (4 4 40) hybrid element are used to discretize the bar. The Lennard Jones pair potential is used to describ e the interactions between atoms. The bar is clamped at one end and a force perturbation is applied to the other end. Figure 4 7 shows the variation of the average stress with time when the stress waves propagat e in the bar. From the curves, it is seen tha t the parallel code is stable. We also compare the current simulation results with those of MD simulation. In Figure 4 7, it is shown that the wave propagation predicted by our method is almost identical to that predicted by MD simulation. Remember that th e degree of freedom here in our FE model is only 1.4% of that of the corresponding MD model. For a successful parallel code with consistency, the computation results should be independent on the number of processors. To test the independence of the simulat ion results on the number of processors, the same computation model is calculated with different number s of processors (1, 4, 8, 16, 32, 64 processors). Figure 4 8 shows that for a varie d number of processors the results are exactly the same. Both Figure 4 7 and 4 8 demonstrate the stability of the parallel code. Scalability Scalability is a measurement that describe s how the parallel code will perform with a varie d number of processors. The scalability test is conducted on the supercomputer Trestle of SDS C ( San Diego Supercomputer Center ). The supercomputer consists of PAGE 67 67 324 computer nodes. Each node contains 32 AMD Magny Cours processors. The clock speed of each processor is 2.4GHz. In this section, two different testing cases are considered: (1) Scaled si ze test (ST): In this test, the number of element s per processor is kept to be 640 With the increasing of number of processors, the size of the whole model is increased proportionally. The total number of elements and the corresponding number of processor s are list ed in Table 4 1. (2) Fixed size test (FT): As shown in Table 4 1, i n this test, the total number of elements is ke pt to be 5120 With an increas ed number of processors, the number of elements assigne d to each processor is reduced. From Figure 4 1, it is known that the total computational time is the sum of the time for steps INPUT, GETMAP, ELIST, EXCHANGE_U, CFORCE, CENTRAL and OUTPUT Among these steps, EXCHANGE_U, CFORCE and CENTRAL a re executed every time step. T he union of these three steps is called SOLVE here. So the total computational time can be calculated using the following equation: (4 1) where N is the number of time steps. The serial fraction is the ratio of computation al time for the serial part to the t otal computational part when using a single processor. For this code, the serial fraction is (4 2) The speedup, which describes how fast the computation would be by using multiple processors, is the ratio of the total computational time by a single processor to the computational time by multiple processors. The expression for speedup, in an ideal situation with no overhead in the parallel part, is PAGE 68 68 (4 3) The efficiency, whi ch describes how efficient the code is when using multiple processors, is defined as the speedup normalized by the number of processors (4 4) In general most of the simulations involve several tens of thousands of time steps. In this work, only one thousand time steps are computed since the speedup and efficiency converged for time steps greater than this number. For all these tests, the computational time for the key steps are listed in Table s 4 2 and 4 3. Usi ng Table 4 2 Ta ble 4 3 and Equation s (4 1) (4 4), the speedup and efficiency for both FT and ST are computed and plot ted in Figure 4 9 and Figure 4 10. It is seen from Figure 4 9 that the speedup for both FT and ST are very close to ideal speedup when the number of pr ocessors is less then 20. As the number of processors continues to increase, the speedup becomes slightly lower than the ideal one. It is worthwhile to mention that the number of elements is another issue which affects the speedup. Fewer elements per proce ssor always leads to higher parallel efficiency. For ST the number of elements assigned to each processor is fixed to be 640. For FT the number of elements assigned to each processor dec r eases with the number of processors used. When 64 processors are us ed in the FT case, there are only 80 elements assigned to each processor. At this time, the communication and other overhead cost for FT is less then that of ST. Therefore, the speedup for FT is higher than for that of ST As shown in Figure 4 10, the effi ciency for both FT and ST drops with an increas ing number of processors. They are both very close to 1 when the number of processors is less then 8. Then as the number of processors increase s the efficiency PAGE 69 69 becomes lower and lower. Fortunately, the rate a t which the efficiency drops becomes lower and lower as the number of processors increase s further Like speedup, the efficiency is also affected by the number of elements assigned to each processor. The efficiency of FT is higher than that of ST. Summary In this chapter, the parallel implementation of the new CG method was reviewed in detail. At first, the four design spaces of the parallel code were discussed Then in the implementation of the code, the overall structure of the parallel code was introduce d T he de tails of each key step of the parallel code were also introduced Finally the performance of the parallel code was tested on the supercomputer Trestle of SDSC ( San Diego Supercomputer Center ). A stress wave propagation example was used for the te st. It was found that the CG simulation results were almost the same as those of the corresponding MD simulation. T he simulation results were independent on the number of processors used. It was also found that the parallel code show ed high efficiency in t he FT and ST tests. As a conclusion, we found that the efficiency of the parallel code is related to the number of elements assigned to each processor. A l arge number of element s per processor always leads to slightly lower efficiency. PAGE 70 70 Table 4 1 Number of elements and processors used in ST and FT ST FT Number of Processors Total number of elements Total number of elements 4 2560 5120 8 5120 5120 16 10240 5120 32 20480 5120 64 40960 5120 Table 4 2 FT results (unit: second) Procs. INP GETM ELIST EXC SOLV OUT 4 2.428 0.179 362.16 0.002 1.647 2.488 8 2.511 0.04 176.05 0.006 0.843 2.493 16 2.538 0.04 88.602 0.006 0.426 2.502 32 2.568 0.05 44.170 0.006 0.220 2.573 64 2.510 0.05 20.071 0.006 0.110 2.4 89 Table 4 3 ST results (unit: second) Procs. INP GETM ELIST EXC SOLV OUT 4 1.226 0.04 163.130 0.002 0.814 1.247 8 2.511 0.04 176.052 0.006 0.843 2.493 16 5.068 0.05 184.924 0.006 0.906 5.062 32 10.021 0.05 188.210 0.007 0.988 10.108 64 19.009 0.05 193.587 0.006 1.076 19.684 PAGE 71 71 Figure 4 1 S tructure of the parallel code Figure 4 2 Procedures for domain decomposition (GETMAP) Figure 4 3 Interfaces between domains PAGE 72 72 Figure 4 4 Communication between neighboring doma ins (EXCHANGE_U) Figure 4 5 Central difference method PAGE 73 73 Figure 4 6 Setup of the model for the force perturbation test Figure 4 7 Vibration of a FCC bar predicted by the two different methods PAGE 74 74 Figure 4 8 Vibration of the FCC bar predicted wi th varies number of processors Figure 4 9 Speedup curves for FT and ST PAGE 75 75 Figure 4 10 Efficiency curves for FT and ST PAGE 76 76 CHAPTER 5 SIMULATIONS OF 2D DY NAMIC BRITTLE FRACTU RE Introduction In this chapter, the proposed coarse grained method (CG AFT) is u sed to simulate the dynamic crack propagation in 2D brittle material. Firstly, the 2D triangular lattice model which has been widely used in the study of 2D brittle fracture is introduced. Then the crack propagation and branching is simulated using both CG and MD method. The comparison of their results are conducted and discussed. At last, as an application of the proposed CG AFT method, five models with different size and loading rate are employed to investigate the interaction of stress waves and crack pr opagation at the atomic scale. As we have mentioned in Chapter 4, the CG AFT method is fundamentally different from MD simulation approaches So the exactly identifying the simulation results by the two different approaches is impractical. But we still can compare their simulation results because both of them are proposed to describe the nature of material at atomic scale. The computational model The first s ignificant MD simulation of crack propagation in a 2D triangular lattice was described by Arhurst and Hoover (1976) in a landmark paper using about 1000 atoms. Large scale MD simulations in parallel computers were performed about 20 years later, for examples, by Holian and Ravelo (1995), Zhou et al. (1996), Abraham et al. (1994, 1996, 1997, 1998, 2000), B uehler et al. (2003), Buehler and Gao (2006). In all of those simulations, the 2D triangular lattice model was again employed, with the computer models containing atoms from thousands to millions of atoms. Harmonic potential s or Lennard Jones potential fun ction w ere employed in most of the simulations PAGE 77 77 to describe the interaction between atoms. A common goal of all these simulations is to model dynamic fracture behavior with general features common to a large class of real physical systems rather than to a p articular material The main objective of this chapter is to demonstrate the application of the new formulation and its finite element implementation. Therefore, in this work, we also use the 2D triangular lattice model in both the FE and the MD simulatio ns. We employ the 2D computer model also out of our consideration of the influence of reflecting wave from the boundaries of the specimen. W ith the available computer power and the capability of our sequential computer code, we wish to generate a computer model that has the dimensions as large as possible so that waves emanated from the crack tip take some time before being reflected back by the boundaries of the specimen to interfere with the forward propagation waves. The well tested and well tuned Lenna rd Jones (LJ) potential used in above mentioned MD simulations is also employed in this work to represent an idealized : (5 1) where r is the separation distanc e between two atoms, is the depth of the potential energy well and is the value of r where the potential become zero. In this work, we employ the normalized parameter system which has been used by Abraham et al. in 1996 and also Holian and Ravelo By u sing this parameter system, all energy quantities are normalized by and all distance quantities are normalized by For the sake of elabration, in the rest of this chapter, we may convert the quantities from this unit system to the physical unit system somewhere, PAGE 78 78 using and Only the nearest neighbors of an atom is T he geometry of the computer mode l is shown in Figure 5 1 The rectangular slab of the 2D triagular lattice has the slab length and width L and W, respectively. A crack with length is cut in the middle of the left boundary. The slab is loaded in model I in the direction perpendicular t o the crack with a constant stain rate. As shown in Figure 5 2, the 2D triangular element introduced in Chapter 3 is used to discretize the model. There are 32,060 triangular 2D elements and 96,180 nodes i n th e FE model shown in the figure. Each element co ntains 55 unit cells. The model corresponds to 1,767,200 atoms (~1.8 million atoms). The degrees of freedom to be solved in the FE model are about 5% of that in the MD model. Simulation details and results Comparison o f CG And MD Simulations In this secti on, two test cases for crack propagation are simulated using both the proposed FE method and MD simulation. Our goal of this section is to compare the results predicted through these two methods. As we have mentioned before, MD simulation has been widely u sed by the pioneers to investigate the brittle crack propagation and branching at atomic scale. Through this comparison, we are trying to show the reliability of the results predicted with the proposed new formula and its finite element implementation. In case 1, the model described in Figure 5 2 is simulated using both of the two methods. The strain loading with the rate of about 10 8 per second (approximately converted from dimensionless unit system to physical unit system) is applied on both PAGE 79 79 the top and the bottom of the slab. The results are shown in Figure 5 3 and Figure 5 4. As predicted in both simulations, at about 1% of the applied strain, the crack begins to open up at the weakest spot, namely the lower tip of the crack or the right corner of the n otch, where the stress is most concentrated. The kink observed in both of the two simulations shows that the crack is unstable at the beginning of its propagation. With the further increasing of the loading, a daughter crack is formed near the main crack. It is seen that the pattern of the crack propagation in the FE simulation does not exactly resemble that in MD simulation. Nevertheless, the FE model does have reproduced the strains for crack initiation and branching similar to that observed in MD simulat ion. Being unable to output the crack velocity history from LAMMPS, we do not output crack speed in this case. However, from the temporal sequence of the snapshot pictures of the propagating crack and the distance it traveled, we see that the crack propaga tion velocities in the FE model are close to that of the MD model. The local stress distribution is also plotted in Figure 5 3. The local stress is calculated using virial stress formula. We understand that the virial stress formula is not appropriate for the description of local stress because it is not consistent with the concept of Cauchy stress in continuum mechanics and does not satisfy the continuum mechanics balance equations (Chen 2006). However since LAMMPS output only local viral stress, for the purpose of comparison we also output the local virial stress in the FE simulation. It is seen from Figure 5 3 that the distribution of the local virial stress are similar in both simulations, but the maximum local stress obtained in the FE simulation is lower than that obtained in the MD simulation. Noticed that MD model contains about 1. 8 millions of atoms and the stress is calculated per atom, while the stress in the FE PAGE 80 80 model is calculated per element and the specimen contain s only 32060 elements. There fore, the MD model has a resolution that is significantly higher than that of the FE model. As a result, the pictures taken from MD simulation have much higher resolution and the maximum stresses are also higher. Figure 5 4 show s the average stress strain curve of the specimen during the processes of tensile loading obtained by the FE and the MD simulations. It is seen that the magnitude of the stress obtained by the FE simulations agrees well with that obtained by MD simulation. However, there is noticeabl e discrepancy between the two stress strain curves, especially between the areas that each curve covers. This implies that the FE model is more brittle than the MD model and underestimates the ductility and the work of fracture of the material. We believe that this considerable discrepancy is caused by the assumption of the linear shape function made in the finite element implementation. In case 2, we use the same geometric model as in case 1. However, the cutoff for interatomic interaction is changed fro m 0.57nm to 0.34nm. The material is thus more brittle than that in case 1. We also reduce the applied strain rate from 10 8 per second in case 1 to 1.810 7 per second in case 2 A temporal sequence of snapshot pictures for the crack dynamics as well as loca l stress distribution by MD and FE are shown in Figure 5 5. Note that we use fixed color bars for both MD and FE results. It is found, this time, that with the more brittle material and slower loading rate the crack start with stable, straight ahead propag ating at the beginning. As the applied strain loading approaching 1.54%, the crack starts to branch. It is worthwhile to meaning that in both of the two cases the crack always starts to propagate at a critical strain about 1%. This PAGE 81 81 means that the loading r ate effect is not so significant to the start of crack propagation as to the start of crack branching. Figure 5 6 shows the local stress distributions in front of the crack tip at the time when the crack starts to propagate. It is seen that the stress dist ributions in these two models are very close. The MD curve consists of 132 points, while there is only 16 points for the FE curve. As a result of the different resolutions in the MD model and the FE model, the stress at the crack tip is lower in the FE mod el. Large Scale FE Simulations We present in this section a series of FE simulations with relatively larger computer models with the same idealized brittle material described in previous section. The goal is to investigate the effect of the interaction of reflected waves with the propagating crack. We simulate the dynamic responses of five specimens that have different size, shape, crack length, and loading rate. The details of the 5 models are given in Table 5 1 Figure 5 7 shows the local stress distribut ions in front of the crack tip for the 5 models measured in the time interval during which the cracks start to propagate. We see that the 5 models predict almost the same critical stress value for crack initiation indepdent of the specimen size and the cra ck lentgh. We also calculate the stress intensity factor, K IC for each model. The results are shown in Table 5 2, from which we find a nearly identical K IC value for all the 5 models. In Figure 5 8, we present the temporal sequence of snapshot pictures of the dynamic process of crack initiation, proagation, and braching for model E, the largest specimen (1.127 m by 0.708 m). We see in Figure 5 8(a) that the stress waves are PAGE 82 82 emitting from the crack tip when the crack starts to propagate, and that the wa ves are undisturbed. However, in Figure 5 8(b), it is seen that the stress waves have been reflected back from the horizontal boundaries and have interfered with the forward propagating crack. As a result, the stress waves are disturbed, and the crack star ts to branch. We see in Figure 5 8(c) that stress waves are more disturbed and the crack has branched into two cracks. In Figure 5 8(d) we observe that each of the two branched cracks further branches into two cracks, which seems to be also resulted from t he interaction of propagating crack with stress waves being reflected back from the horizontal boundaries. Two newly branched cracks soon arrest, leaving two main cracks continue to propagate. This is why we see kinks in Figure 5 8(e) and (f), and the two main bra nches propagate to the right boundary of the specimen in a zigzag pattern. We plot the time history of the crack speed in Figure 5 9 for the period of crack initiation and first branching, in which is the Rayleigh wave s peed. It is seen that the first sharp increase in crack speed occurs when the reflected waves meet with the forward propagating stress waves at the strain around 0.58% (cf. Figure 5 8(b)). Then, we observe a sharp drop of the crack speed after the branchi ng (cf. Figure 5 9). Point A and B correspond to the crack initiation and the first crack branching We see before and after branching events the crack speed oscillates at about 0.4 Rayleigh speed From Figure 5 8 and Figure 5 9 we find that the sudden increase of crack speed is caused by the interaction of forward stress waves with reflected stress waves, which in turn cause the crack to branch, namely the dynamic instability of brittle fracture. Figure 5 10 present s the snapshot picture s taken at appr oximately the time of the first c rack branching for all the five models. The left are zoomed in pictures at the crack PAGE 83 83 tips. Again, it is seen that cracks, in all the five models, branch as a direct result of the interaction between the reflected stress wav es and the forward propagating cracks Table 5 3 shows the time for different waves to emit from and travel back to the crack tip and also the time interval (the output time interval in the FE simulation) during which the crack branching occurs. Here, and stands for the time for longitudinal waves, shear waves, and Rayleigh waves to emit from and travel back to crack tip, respectively. It is found that in all of the five models the longitudinal waves always arrive at the crack tip slightly before the crack branching. The shear waves and Rayleigh waves arrive at the crack tip after the crack branching, and in all the five models it is the longitudinal waves that causes the crack speed to rapid increase and subsequently the cracks to branch. Summary We have presented in this chapter a series of 2D simulation of the dynamic fracture. Using the 2D triangular linear finite element introduced in Section 3.2 a 2D computational model w ith about 5% degrees of freedom of the corresponding MD model was constructed. By comparing the coarse grained and MD simulation results, we have showed that the coarse grained method based on the atomistic formulation can predict the dynamic evolution of the crack tip stress, crack propagation speed, crack path and branching, qualitatively and quantitatively similar with atomic level molecular dynamics simulations. Form the larger scale model (~7 millions of atoms), we have found that the interaction of st ress waves emitting from the propagating crack tip and the stress waves reflected back form the boundaries has a strongly influence on the crack speed, crack path and crack branching, PAGE 84 84 Besides the given temperature limitation, the linear shape function inte rpolation used here is also a factor that affects the accuracy of the computations. As we know, at the crack tip the displacement field always has very high gradient. Thus, some more sophisticate shape functions are needed to improve the accuracy of the co mputation. PAGE 85 85 Table 5 1. Geometric parameters and loading rates for the six models Models A B C D E Specimen l ength (Lx) (nm) 1127 1127 1127 566 1127 Specimen w idth (Ly) (nm) 354 354 354 354 708 Crack length (a) (nm) 56 225 450 225 450 Number of nodes 38 4,480 384,480 384,480 192,480 768,960 Number of elements 128,160 128,160 128,160 64,160 256,320 Number of atoms 7,048,800 7,048,800 7,048,800 3,528,800 14,097,600 Applied strain rate (1/s) 4.610 6 4.610 6 4.610 6 9.210 6 4.610 6 Table 5 2. K IC for mo dels A~E Models A B C D E K IC ( ) 1.77 1.71 1.71 1.73 1.65 Table 5 3. Time interval during which crack branching occurs and time of arrival of reflected stress waves Models A B C D E Time interval for branching (ps) 80~ 120 16 0~ 2 0 0 160~ 2 0 0 160~20 0 320~ 3 6 0 Time of arrival of Longitudinal wave (ps) 85 243 464 249 828 (ps) 176 176 176 176 352 Time of arrival of Shear wave (ps) 147 421 803 431 1433 (ps) 305 305 305 305 609 Time of arrival of Rayleigh wave (ps) 158 453 863 463 1541 (ps) 328 328 328 328 655 PAGE 86 86 Figure 5 1 Simulation Model and Boundary Conditions Figure 5 2 The finite element mesh of the computer model PAGE 87 87 Figure 5 3 Comparison between FE and MD simulations results (deformation and stress distribution) for case 1 Left: MD simulation results (~1.8 million atoms), Right: FE simulation results (~32K elements) PAGE 88 88 Figure 5 4 Average stress vs. strain curves Figure 5 5. Comparison of FE and MD simulations results (deformation and stress distribution) for case 2 Left: MD simulation results (~1.8 million atoms), Right: FE simulation results (~32K elements) PAGE 89 89 Figure 5 5 Continue d PAGE 90 90 Figure 5 6 Stress distribution in front of crack tip when crack start to propagate Figure 5 7 Stress distribution in front of crack tip when crack start to propagate PAGE 91 91 Figure 5 8 FE results of dynamic crack propagation. The pictures a re taken at strain (a) 0.49%, (b) 0.58%, (c) 0.75%, (d) 0.97%, (e) 1.07%, (f) 1.18% (a) ( b ) ( c ) PAGE 92 92 Figure 5 8 Continue d ( d ) ( e ) ( f ) PAGE 93 93 Figure 5 9 Crack speed measured for model E during the period of crack initiation and first branching Figure 5 10 Stress waves at the time of crack branching in models A E ( Left: zoomed in the crack tip ). Model A Model B Model C PAGE 94 94 Figure 5 10 Continue d Model D Model E PAGE 95 95 CHAPTER 6 SIMULATIONS OF 3D DY NAMIC BRITTLE FRACTU RE Introduction Three d imensional dynamic fractures are much more complex than two dimensional cases not only because of the addition of one dimension but also because of the new phenomen a that would be involved with the addition of the dimension. On the other hand, the analyti cal solution for 3D dynamic fracture problem s is always too difficult to derive. Thus a numerical simulation tool for 3D dynamic fracture is necessary. In the prev ious chapter, we used the new CG method to study 2D dynamic fracture problems. Because of the differences between 2D and 3D problems, we still need to evaluate the performance of the new CG method in simulating 3D dynamic fracture. Firstly, we compare the re sults of CG and MD simulations to show that the results computed by the CG method are reas onable physic ally It is worthwhile to mention that there is a difference between the treatments of the temperature field in the CG and the MD simulations. In the CG simulation, temperature is an independent variable which can be determined by solving the conservation of energy equation. In this work, as a special case, temperature is kept to be 0 K throughout the CG simulation. While in the MD simulation under the NVE ensemble, temperature is coupled with the velocity of each atom in the system. Even if we start with a homogeneous 0 K temperature, the temperature could be very high and inhomogeneous ly distributed after the simulation. The NVT or NPT ensemble is always employed for the sake of controlling the average temperature. However, neither the NVT nor NPT algorithm is suitable for the comparison here due to the numerical dumping that they would bring. This dumping can significantly change the wave speed in the model. Therefore i n this chapter, we chose PAGE 96 96 the NVE ensemble in the MD simulation. When the system temperature in the MD simulation is low, the two simulation results are comparable. Next, a series of CG simulation s including 3D crack branching, impact test s and 3D crack surface evolution, are performed Through these simulations, we want to sh ow that the new CG method can predict the arbitrary dynamic crack propagation in crystals. In the last example, 512 processors is used for a system composed of 0.11 million elements. With this example, the performance of the parallel code in handling large scale simulations is also demonstrated. Comparison of CG and MD simulations Crack propagation is associated with the loss of material stability which ushers in a host of numerical challenges. In this section, the 3D crack propagation in an FCC single crys tal film is simulated using both CG and MD methods. Through the comp arison of their results, we want to show the reliability of the CG method in pred icting 3D dynamic fracture. In this example the robustness of the CG method is demonstrated through the si mulation of 3D dynamic fracture As shown in Figure 6 1 the computational model is a notched thin film which is discretize d into 1590 hybrid elements T he thin film is fixed at the bottom and a constant velocity 20m/s, is applied to the top surface in th e direction shown by the arrow. For the purpose of comparison, a n MD model with the same size, shape and boundary conditions is computed using LAMMPS. Since each hybrid element used here is equivalent to 1000 atoms, the MD model contains 1,590,000 atoms. F igure 6 2 (a) and (b) show local stress ( ) distribution s for the two models PAGE 97 97 when the average strain is 1.9%. At this stage, there is no crack propagation and the material stability is still maintained Both of the two pictures ( Figure 6 2 (a) and (b) ) show that there are two stress concentration regions. One is at the crack tip ; the other is at the top right corner of the model. The comparison of these two pictures shows that, before crack propagation, t he CG model predicts quite similar local material behavior to the corresponding MD model. When the crack start s to propagate, because of the stress concentration, the crack path will bend up to the top right corner of the model As mentioned in section 3, the crack propagation always prefers the weak surfaces in crystals. The results shown in Figure 6 2 ( c ) and ( d ) confirm this point Both the CG result ( Figure 6 2 ( c ) ) and the MD result ( Figure 6 2 ( d ) ) show that the crack path bends up to the top right re gion of the model with an angle of 60 o This 60 o angle corresponds to one of the slip planes of the material. Note that the differences between the CG result and the MD result are caused by the different treatment of temperature. As mentioned before, the s ystem temperature is kept to 0 K throughout the CG simulation. Before crack propagation ( 1.9%) the temperature of the MD model is less than 10 K It is reasonable that the MD result is similar to the CG result at such a low tempe rature. After crack propagation, the temperature of the MD model increase s suddenly. At the stage of =2.9%, the temperature of the MD model is over 100k. So the CG result is different from the MD result as shown in Figure 6 2 ( c ) and ( d ) Figure 6 3 shows the average stress strain curves calculated by the CG and the MD simulations It is found that the two curves are very close to each other before they PAGE 98 98 reach their respective peaks corresponding to the beginning of crack propagatio n. T he staircases of these two curves are caused by stress wave propagation in the specimen After crack propagation, the temperature of the MD model is much higher than that of the CG model. Thus the fracture in the CG model is more brittle because of the lower temperature. This is the reason why the CG curve drops earlier than the MD curve It is worthwhile to mention that the crack propagation, especially when it becomes unstable, is ver y sensitive to temperature Since the CG and the MD simulations are p erformed under different thermodynamic conditions the crack paths predicted by the two m ethods are not exactly the same To obtain close r results, temperature should not be assumed to be a constant in the CG simulation. This will lead to the addition of a n unknown to the computation. To solve this unknown, the conservation of energy equation need s to be solved in combination with the balance of linear momentum equations. Application of the CG method in 3D brittle fracture 3D Crack Branching Under Mixed Mod e Loading In this section, the crack branching in an FCC single crystal is simulated The computational model is shown in Figure 6 4 (a) (c) The thin film is discretized into 8,944 hybrid elements which is equivalent to 8.944 milion atoms The dimensio ns of the thin film are 325.4nm 108.5nm 1.8nm. Figure 6 4 (a) shows the front view of the model. A notch with the length of 130.2nm ( 40% of the length of the thin film) is cut in the middle of the film. is the constant distribute d force applied to the top and the bottom surfaces Figure 6 4 (b) is the side view of the model. Figure 6 4 (c) is the magnified view of (b). From Figure 6 4 (b) and (c), it is seen that the loading will cause mixed mode fracture in the specimen. The CG s imulation is performed using 16 PAGE 99 99 processors. For this CG model, it is assumed that the crack propagation and branching always prefer the slip planes. This assumption is reasonable according to Dr. Gumbsch and his co They found that the cr ack branching in FCC single crystal s always follows the path of dislocations emitted from the crack tip (Gumbsch et al. 1997). They also explained that the dislocation nucleation distorts the arrangement of the atoms near the slip plane at the crack tip an d creates a weak path for the crack to follow. In my work, a n MD simulation is performed to demonstrate the validity of this conclusion for a mixed mode loading case. Because of the high cost of the MD simulation, t he model is reduced to about one eighth o f the CG model. Each dimension of the MD model is about half of that of the CG model s o there are about 1million atoms in the MD m ode l The MD simulation ( Figure 6 5 ) shows that the cr ack start s to propagate stably at the beginning Then it branches into the directions of +60 o and 60 o relative to the original path. These directions correspond to the preferred directions of dislocations. This result shows that the conclusion of Dr. Gumbsch and his co workers is also valid for mixed mode fracture problems As shown in Figure 6 6 the simulation result of the larger CG model is similar to the MD result shown in Figure 6 5 Although the size of the CG model and the MD model are different, the similar crack propagation pattern and the crack branching angle have revealed the nature of the crack branching in this problem. T he dislocations emitted from the crack tip determine the paths for crack branching in this material. Impact Fracture o f a Plate In this case, a more complex fracture problem is computed using t he CG method T he applicability of the method to large deformations, fragmentation and multiple crack PAGE 100 100 problems is conducted T he computational model for this impact test is shown in Figure 6 7 A rigid ball with a constant velocity V hit s on the FCC single crystal plate The FCC plate contains 40 000 hybrid elements which are equivalent to 40,000,000 atoms. This plate is a very thin film with the dimension s of 511nm 511nm 1.8nm The rigid ball has a radius of 50nm. It moves towards the plate with a constan t velocity of 100m/s. The simulations in this test were conducted with 64 processors. The snapshots of this impact simulation are shown in Figure 6 8 The contour represents the displacement field DZ in the deformed body. DZ indicates the displacement norm al to the plate surface and opposite to the moving direction of the rigid ball. As shown in Figure 6 8 when the rigid ball hit the plate, the plate fractured with six small cracks Then the six cracks propagate in six directions with almost the same crack speed The angle bet ween each two neighboring crack branches is nearly 60 o The bottom right picture in Figure 6 8 shows the oblique view of the deformed specimen at t=200ps. Because of the large size of the CG model, the corresponding MD simulation is t oo computationally expensive to be performed Therefore i n this work, a similar impact simulation on a much smaller model is conducted with MD. The MD model has the same thickness as the CG model. But the dimensions in the other two dimensions are only on e tenth those of the CG model. Correspondingly the size and the velocity of the rigid ball are also reduce d to one tenth of the CG model. Although the size and the shape for these two models are different, their simulation results reveal the nature of the problem. From Figure 6 9 although the number of cracks is reduced to four the fracture pattern of the plate is still similar to that in the CG result The four cracks propagate PAGE 101 101 symmetrically and with almost the same velocity. The reason for the reductio n in the number o f cracks in the MD simulation is that the ratio of the thickness to the in plane dimensions is ten times larger than that of the CG model. Since thicker of plate s dissipate more energy from an impact than do thinner plates only four crack s are sufficient to dissipate the energy brought by the impact of the rigid ball. 3D Crack Surface Evolution The objective of this example is two fold. Firstly, we want to simulate 3D crack propagation in a single crystal block using the CG method Secondl y, we want to test the performance of the parallel code for large scale simulation. The computational model is s hown in Figure 6 1 0 A notched FCC rhombohedra l block is fixed on the bottom and a constant surface loading is app lied to the top of the model. T he applied loading is set to be much larger than that for stable crack propagation. This computational model is composed of 1,545,348 nodes and 110,382 hybrid elements which correspond to over 0.11 billion atoms. 512 process ors are employed for the computation to test the performance of the parallel code in handling a large number of processors Figure 6 1 1 shows some snapshots of this complex c rack surface propagation process The contour in the pictures indicates the displa cement field DY in the specimen. DY is the displacement along the direction of the applied force and normal to the top surface of the model. The snapshots show the progress of the 3D crack surface evolution under the effects of both stress waves and the cr ystal structure of the material. It is seen that the crack surface is smooth at the beginning of the crack propagation. Then because of the instability of the crack propagation the crack surface becomes rougher and rougher. From the snapshots shown in Fig ure 6 1 1 it is also found that the PAGE 102 102 crack surface is smoother when close to the boundary of the specimen and rougher in the middle of the specimen This different crack propagation pattern indicates the complex stress field within the 3D model. Summary We presented in this chapter a series of 3D dynamic fracture simulations using the hybrid element introduced in C hapter 3 As we know, losing part of the atomic information is inevitable for coarse graining. Then the question is how much confidence we would have when ignoring such atomic information. As we have mentioned in the first chapter, the fracture behavior of a single crystal material is close ly related to its crystal structure. So it is necessary to maintain as much atomic information as we can i n t he CG model to obtain the most accurate results In this chapter, it was shown that the utilization of the hybrid element is a possible way to achieve this goal. By comparing the CG simulation results with the MD simulation results, we have shown that the new coarse grained method can predict the 3D crack tip stress, stable crack propagation, and unstable crack propagation both qualitatively and quantitatively similar with the atomic level MD simulations It is worthwhile to mention that the difference bet ween coarse grained and MD results comes from the different thermodynamic processes involved in the two different simulations To obtain close r results for the two methods the conservation of energy equation needs to be solved along with the balance of li near momentum equations. As applications of the CG method, 3D crack branching under mixed mode loading, 3D impact test, and 3D unstable crack surface evolution were simulated to show the robustness of the coarse grained method in dealing with large scale 3 D dynamic fracture. In the last case, 3D crack surface evolution under M ode I loading, 512 PAGE 103 103 processors were employed to compute a model with over 0.1 million elements which is equivalent to over 0.1 billion atoms in the corresponded MD model. Through this case, the ability of the parallel code in handling a large scale model with a large number of processors was demonstrated. PAGE 104 104 Figure 6 1 Notched film for the first test case (a) CG (b) MD (c) CG (d) MD Figure 6 2 Comparison of CG and MD results for the first case yy = 1.9% yy = 2.9% PAGE 105 105 Figure 6 3 Comparison of the average stress time curves Figure 6 4 Computation model for 3D crack branching (b) (a) (c) PAGE 106 106 Figure 6 5 MD simulation of c rack branching under mixed mode loading Figure 6 6 CG sim ulation of c rack branching under mixed mode loading. PAGE 107 107 Figure 6 7 Computation model for 3D impact test Rigid ball Thin plate composed of 40,000 hybrid elements PAGE 108 108 Figure 6 8 Fracture of the plate under the impact of a rigid ball Figure 6 9 Impact test for a smaller plate by MD simulation PAGE 109 109 Figure 6 1 0 Computation model for the bulk material under M ode I loading Figure 6 1 1 Crack surface evolution in bulk material PAGE 110 110 Figure 6 1 1 Continue d PAGE 111 111 Figure 6 1 1 Continue d PAGE 112 112 Figure 6 1 1 Continue d PAGE 113 113 CHAPTER 7 SUMMARY AND FUTURE W ORK Summary and conclusion The simulation of dynamic fracture is one of the most challenging topics in solid mechanics. On one hand, the existence of a crack or dislocation always brings discontinuities to the problem. On the other hand, fracture, in nature, is a multisca le dynamic process which involves the break ing of material at different length scales rang ing from the macroscale to nanoscale. Molecular dynamics is a numerical approach that successfully meets the two challenges list ed above. However, this approach can only handle a model with a very small size because of the large time and memory cost. Multiscale modeling is a promising solution, but currently there is no successful model which deal s with large scale dynamic fracture. In this work, a new coarse grained method which incorporates almost all of the critical atomic information embedded in to large elements, was presented. It is well known that losing part of the atomic information is inevitable for coarse graining. Then the question is how much confidence we have when ignoring such atomic information. As we have mentioned in the first chapter, the fracture behavior of a single crystal material is close ly related to its crystal structure. Thus not only the interatomic potential but also the key feature of the crystal structure should be embedded in the coarse grained model. Serving as the constitutive relation in AFT, the interatomic potential is naturally embedded in the coarse grained model b ut how can the key features of the crystal structure in the coarse grained model be maintained ? In this work, the utilization of carefully designed finite elements provides a prospective way to achieve this goal. For a 2D triangular lattice model, a 2D triangular finite element is designed to capture PAGE 114 114 the main feature of t he lattice structure. For a 3D single crystal with a rhombohedral primitive cell, a hybrid finite element aimed to deal with the crack or dislocation propagation is designed. To improve the efficiency and extend the size of the computation model, a paralle l code is implemented using FORTRAN combined with MPI (Message Passing Interface). Domain decomposition strategy is used to divide the entire model into several sub models. Each sub model is assigned to one processor a nd the communication between neighbor ing processors is performed every time step to share part of their results. The performance of the parallel code is tested on the supercomputer Trestle of SDSC (San Diego Supercomputer Center). The results show that the parallel code can handle large scale computation with high efficiency. To validate the coarse grained method, the comparison of the coarse grained results and the correspond ing MD results are conducted. The results show that the new method can predict the crack tip stress, dynamic crack prop agation and even crack branching, all of which are qualitatively and quantitatively similar with atomic level molecular dynamics simulations in both 2D and 3D cases. It is worthwhile to mention that the differences of the results do not solely come from th e coarse graining. The temperature treatment is also an important source of the differences. To approximate the reality more accurately, the complex temperature field around the crack tip needs to be considered. As applications of the developed coarse gr ained method, crack branching is observed in both 2D and 3D cases. In the 2D case, it is found that the interaction of stress waves emitting from the propagating crack tip and the stress waves reflected PAGE 115 115 back from the boundaries ha ve a strong influence on t he crack speed, crack path and crack branching. In 3D case, it is found that crack branching is close ly related to the dislocations nucleated near the crack tip. These dislocations create two weak paths for the crack to follow. In the 3D simulation of an impact test on a thin plate, the applicability of the method to large deformation, fragmentation and the propagation of multiple unstable cracks is demonstrated. In the case of 3D unstable crack surface evolution, the ability of the parallel code in handl ing large scale computation with a large number of processors is shown. Future work In this work, we only consider the special case where the temperature field in the whole model is given in advance of the computation. Sometimes this is a very crucial cond ition. For example, the temperature distribution at the crack tip is always very complex with a high gradient due to the local heat generation during dynamic crack propagation. This means that the dynamic crack propagation phenomenon is always a thermal/m echanical coupling problem. So the conservation of energy equation should be solved in combination with the balance equation of linear momentum. In AFT the velocity of each atom is divided into two parts: the stream velocity and the instantaneous velocity The difference is the velocity which contributes to the temperature of the system. In the presented coarse grained method, the stream velocity field within an element is approximated by linear shape functions. This approximation is perfect for homogeneou s deformation. However, when the deformation gradient is large, such as around the crack tip, this approximation may be too simple if using a very large element. To overcome this element size limitation, higher order shape functions or more PAGE 116 116 sophisticated s hape functions should be employed. The presented methodology is a promising tool for the investigation of the mechanism of dynamic fracture of crystalline material s The scope of the application should not be confined to some simple crystal such as FCC or BCC single crystal s In our future work, this methodology could be used in coarse graining the atomistic dynamics of many complex materials such as: semiconductor materials, biomineral materials, or even crystalline polymer materials. PAGE 117 117 LIST OF REFERENCE S Abraham F.F., 1996. Dynamics of Brittle Fracture with Variable Elasticity Physical Review Letters 77, 869. Abraha m F.F., Brodbeck D., Rafey R.A., Rudge W.E., 1994 Instability dynamics of fracture: A computer simulation investigation. Physical Review Letters 73, 272 275. Abraham, F.F., Brodbeck, D., Rudge, W.E., Broughton, J.Q., Schneider, D., Land, B., Lifka, D., Gerner, J., Rosenkrantz, M., Skovira, J. Gao, H. 1998. Ab initio dynamics of rapid fracture. Modelling and Simulation in Materials Science and Engineering 6, 639 670. Abraham, F.F., Gao, H., 2000. How Fast Can Cracks Propagate? Physic al Review Letters 84 3113 Anderson, T.L., 1991. Fracture Mechanics, Fundamentals and Applications. CRC. Broughton, J.Q., Abraham, F.F., Bernstein, N., Kaxiras, E., 1999. Concurrent coupling of length scales: methodology and application. Physical revie w B 60, 2391. Buehler, M.J., 2010. Strength in numbers. Nature Nanotechnology 5, 172 174 Buehler, M.J., Abraham, F.F., Gao, H., 2003. Hyperelasticity governs dynamic fracture at a critical length scale. Nature 426, 141 146. Buehler, M.J., Gao, H., 20 06. Dynamical fracture instabilities due to local hyperelasticity at crack tips. Nature 439,307 310. Camacho, G.T., Ortiz, M., 1996. Computational modeling of impact damage in brittle materials. International Journal of Solids and Structures 33, 2899 93 8. Chen, Y., 2006 Local stress and heat flux in atomistic systems involving three body forces. Journal of Chemical Physics 124, 054113. Chen, Y., 2009. Reformulation of microscopic balance equations for multiscale materials modeling. Journal of Chemic al Physics 130, 134706 Chen, Y Lee, J.D., 2005. Atomistic Formulation of A Multiscale Theory for Nano/Micro Physics. Philosophical Magazine 85 4095. Chen, Y., Lee, J.D., Lei, Y., Xiong, L., 2006. A Multiscale Field Theory: Nano/Micro Materials, G. C. Sih (Eds), Book chapter (invited contribution) on Multiscaling in Molecular and Continuum Mechanics Springer, 23. PAGE 118 118 Coffman, V.R., Sethna, J.P., Heber, G., Liu, M., Ingraffea, A., Bailey, N.P., Barker, E.I., 2008. A comparison of finite element and atom istic modeling of fracture. Modelling and Simulation in Materials Science and Engineering 16, 065008. Cox, B.N., Gao, H., Gross, D., Rittel, D., 2005. Modern topics and challenges in dynamic fracture. Journal of Mechanics and Physics of Solids 53, 565 5 96. Deng, Q., Xiong, L., Chen, Y., 2010. Coarse graining atomistic dynamics of brittle fracture by finite element method. International Journal of Plasticity 26, 1402 1414. Deng, Q., Chen, Y., Lee, J., 2009. An Investigation of the Microscopic Mechanism of Fracture and Healing Processes in Cortical Bone. International Journal of Damage Mechanics 18, 491 502. Eringen, A.C. 1972a. Linear theory of nonlocal elasticity and dispersion of plane waves. International Journal of Engineering Science 10, 425 43 5. Eringen, A.C. 1972 b Nonlocal polar elastic continua. International Journal of Engineering Science 10, 1 16. Eringen, A.C. Edelen, D.G.B., 1972. On nonlocal elasticity. International Journal of Engineering Science 10, 233 248. Eringen, A.C., Suh ubi, E.S., 1964. Nonlinear theory of simple micro elastic solids I. International Journal of Engineering Science 2, 189 203. Espanol, P., 2004. Statistical mechanics of coarse graining. Lecture Notes in Physics 640, 69. Flores, K.M., Dauskardt, R.H., 1 999. Enhanced Toughness Due to Stable Crack Tip Damage Zones in Bulk Metallic Glass. Scripta Materialia 41, 937 943. Flynn, M.J., 1972. Some computer organizations and their effectiveness. IEEE Transactions on Computers C 21(9). Freund, L.B., 1974. Cr ack propagation in an elastic solid subject to general loading. IV. Obliquely incident stress pulse. Journal of Mechanics and Physics of Solids 22, 137 146. Freund, L.B., 1998. Dynamic Fracture Mechanics 2 nd edn. Cambridge Univ. Press, Cambridge, UK. Fr ies, T.P., Belytschko, T., 2006. The intrinsic XFEM: a method for arbitrary discontinuities without additional unknowns. International Journal for Numerical Methods in Engineering 68, 1358 1385. PAGE 119 119 Gao, H., 1994. Surface roughening and branching instabiliti es in dynamic fracture. Journal of Mechanics and Physics of Solids 41, 457 486. Ghoniem, N.M., Busso, E.P., Kioussis, N., Huang, H., 2003. Multiscale modeling of nanomechanics and micromechanics: an overview. Phylosophy Magzine 83, 3475 3528. Gibbs, J. W., 1902. Elementary principals in statistic mechanics. Yale University Press. Gumbsch, P., Zhou, S.J., Holian, B.L., 1997. Molecular dynamics investigation of dynamic crack stability. Physical Review B 55, 3445 3455. Hardy, R.J., 1982. Formulas for det ermining local properties in molecular dynamics simulation: Shock waves. J. Chem. Phys. 76, 622 628. Holian B.L., Ravelo R., 1996. Fracture simulations using large scale molecular dynamics. Physical Review B 5, 11275. Hone, T.D., Izvekov, S., Voth, G.A., 2005. Fast centroid molecular dynamics: A force matching approach for the predetermination of the effective centroid forces. J. Chem. Phys. 122, 54105. Hoover, W.G., 1986. Molecular Dynamics. Springer, Berlin. Hoover, W.G., 1991. Computational Statistical Mechanics. Elsiver, Amsterdam. Huynh, D.B.P., Belytschko, T., 2009. The extended finite element method for fractu re in composite materials. International Journal for Numerical Methods in Engineering 77, 214 239. Irving, J.H., Kirkwood, J.G., 1950. The statistical mechanics theory of transport process: IV. The equations of hydrodynamics. J. Chem. Phys. 18, 817. Ir ving, J.H., Zwanzig, R.W., 1951. The statistical mechanics theory of transport process: V. The quantum dynamics. J. Chem. Phys. 19, 1173. Izvekov, S., Parrinello M., Burnham C.J., Voth G.A., 2004. Effective force fields for condensed phase systems from a b initio molecular dynamics simulation: A new method for force matching J. Chem. Phys. 120, 10896. Izvekov, S., Voth, G.A., 2005. Effective Force Field for Liquid Hydrogen Fluoride from Ab Initio Molecular Dynamics Simulation Using the Force Matching Me thod. J. Chem. Phys. B 109, 6573. Kirkwood, J.G., 1946. The statistical mechanics theory of transport process: I. General theory. J. Chem. Phys. 14, 180. PAGE 120 120 Kirkwood, J.G., 1947. The statistical mechanics theory of transport process: II. Transport in Gass es. J. Chem. Phys. 15, 72. Kirkwood, J.G., Buff F.P., Green, M.S., 1949. The statistical mechanics theory of transport process: III. The coefficient of shear and bulk viscosity of liquids. J. Chem. Phys. 17, 988. Lyubartsev, A.P., Laaksonen, A., 1995. Calculation of effective interaction potentials from radial distribution functions: A reverse Monte Carlo approach. Phys. Rev. E 52, 3730. Lyubartsev, A.P., Laaksonen, A., 1996. Effective potentials for ion DNA interactions. J. Chem. Phys. 11, 207. Mat tson, T.G., Sanders, B.A., and Massingill, B.L., 2005. Patterns for Parallel Programming. Pearson Education. Miller, R., Ortiz, M., Phillips, R., Shenoy V.B., Tadmor, E.B., 1998. Quasicontinuum models of fracture and plasticity. Engineering Fracture Mecha nics 61, 427 444. Needleman, A., 1987. A continuum model for void nucleation by inclusion deboning. ASME Journal of Applied Mechanics 54, 525 531. Park K., Paulino G.H., Roesler J.R., 2009. A unified potential based cohesive model of mixed mode fractur e. Journal of Mechanics and Physics of Solids 57, 891 908. Pashley, D.W., 1960. A Study of the Deformation and Fracture of Single Crystal Gold Films of High Strength inside an Electron Microscope. Proceedings of the Royal Society of London. Series A, Mat hematical and Physical Sciences 255,218 231. Pons, A.J., Karma, A., 2010. Helical crack front instability in mixed mode fracture. Nature 464, 85 89. Remmers, J.J.C., Borst, R., Needleman, A., 2003. A cohesive segments method for the simulation of crack growth. Computational Mechanics 31, 69 77. Ross, J., Kirkwood, J.G., 1954. The statistical mechanics theory of transport process: VI. Quantum theory of transport in gasses. J. Chem. Phys. 22, 1094. Schultz, R.A., Jensen, M.C., Bradt, R.C., 1994. Singl e crystal cleavage of brittle materials. International Journal of Fracture 65, 291 312. Silling, S.A, 2000. Reformulation of Elasticity Theory for Discontinuities and Long Range Forces. Journal of the Mechanics and Physics of Solids 48, 175 209. PAGE 121 121 Sillin g, S.A, Ebrahim, A., 2005. A Meshfree Method Based on the Peridynamic Model of Solid Mechanics. Computers and Structures 83, 1526 1535. Soper, A.K., 1996. Empirical Monte Carlo simulation of fluid structure. Chem. Phys. 202, 295. Soper, A.K., 1997. The quest for the structure of water and aqueous solutions. J. Phys.: Condens. Matter 9, 2717. Tadmor, E.B., Ortiz, M., Phillips, R., 1996. Quasicontinuum analysis of defects in solids. Philosophy Magzine 73, 1529 1563. Xiao, S.P., Belytschko, T., 2004. A bridging domain method for coupling continua with molecular dynamics. Computer Methods in Applied Mechanics and Engineering. 193, 1645 1669. Xiong, L., Chen, Y., 2009a. Coarse grained simulations of single crystal silicon. Modelling and Simulation in Mat erials Science and Engineering 17, 035002:1 035002:17. Xiong, L., Chen, Y., 2009b. Multiscale modeling and simulation of single crystal MgO through an atomistic field theory. International Journal of Solids and Structures 46, 1448 1455. Xiong, L., Chen Y., Lee, J., 2008. Simulation of dislocation nucleation and motion in single crystal magnesium oxide by a field theory. Computational Materials Science 42, 168. Xiong, L, Tucker, G., McDowell, D.L., Chen, Y., 2011. Coarse grained atomistic simulation o f dislocations. Journal of the Mechanics and Physics of Solids 59, 160 177 Yamakov, V., Saether, E., Phillips, D.R., Glaessgen, E.H., 2006. Molecular dynamics simulation based cohesive zone representation of intergranular fracture processes in aluminum. Journal of the Mechanics and Physics of Solids 54, 1899 1928. Yoffe, E.H., 1951. The moving Griffith crack. Philosophy Magzine 42, 739 750. Zhou S.J., Lomdahl P.S., Thomson R., Holian B.L., 1996. Dynamic Crack Processes via Molecular Dynamics. Physical Review Letters 76, 2318. Zwanzig, R.W., Kirkwood, J.G., Stripp, K.F., Oppenheim, I., 1953. The statistical mechanics theory of transport process: VI. A calculation of the coefficients of shear and bulk viscosity liquids. J. Chem. Phys. 21, 2050. PAGE 122 122 Zwanzig, R.W., Kirkwood, J.G., Oppenheim, I., Alder, B.L., 1954. The statistical mechanics theory of transport process: VI. The coefficient of thermo conductivity of monatomic liquids. J. Chem. Phys. 22, 783. PAGE 123 123 BIOGRAPHICAL SKETCH Qian Deng w as born in Jingzhou Hubei P.R. China in 1978. After spending his early years in Jingzhou Qian moved to Dalian, P.R. China in 1996 to attend the Dalian University of Technology (DUT) He received his B.S. in Engineering Mechanics from DUT in 2000. Then h e moved to Wuhan, Hubei, P.R.China to pursue a M.S. in Huazhong University of Science and Technology (HUST). He received his M.S. in Engineering with a major in Solid Mechanics from HUST in 2006. At the same year oup to pursue a Ph.D. in solid mechanics H e firstly spend one year to study with the mentorship of Dr. James Lee at the George Washington University (GWU) Then, following the plan, he joined the Atomistic and Multiscale Mechanics lab at the department of Mechanical and Aerospace Engineering, University of Florida. He received his Ph.D. from the University of Florida in the fall of 2011. His current research focuses on developing a multiscale p arallel simulation tool for applications in computational material science 