

small (250x250 max)
medium (500x500 max)
Large
Extra Large
large ( > 500x500)
Full Resolution


MULTISCALE SIMULATION FROM ATOMISTIC TO CONTINUUM  COUPLING MOLECULAR DYNAMICS (MD) WITH MATERIAL POINT METHOD (MPM) By NITIN P. DAPHALAPURKAR Bachelor of Engineering University of Pune Pune, Maharashtra, India 2002 Submitted to the faculty of the Graduate College of the Oklahoma State University in partial fulfillment of the requirements for the Degree of MASTER OF SCIENCE December 2004 ii MULTISCALE SIMULATION FROM ATOMISTIC TO CONTINUUM  COUPLING MOLECULAR DYNAMICS (MD) WITH MATERIAL POINT METHOD (MPM) Thesis Approved: Dr. Hongbing Lu Dr. Ranga Komanduri Dr. Samit Roy Dr. A. Gordon Emslie Advisor Dean of the Graduate College iii ACKNOWLEDGEMENTS First, I would like to express my sincere gratitude and appreciation to Professor, Dr. Hongbing Lu for his guidance, support, and encouragement during the entire course of this study. Working under his advisement was really a valuable experience for me. I would also like to express my appreciation to Dr. S. Roy and Dr. R. Komanduri for their constant guidance and encouragement throughout this study. Thanks are due to Dr. K. A. Salam for reviewing my thesis work. I would like to express my appreciation to Dr. Bo Wang for his invaluable guidance and direct contribution towards the completion of this thesis. Without their great friendship, understanding, and encouragement, it would be hard to believe that I could come this far. My master studies under them have been a wonderful and unforgettable period. The work was supported by a grant from the Air Force Office of Scientific Research (AFOSR) through a DEPSCoR grant (No. F496200310281). The author thanks Dr. Craig S. Hartley, Program Manager for Metallic Materials Program at AFOSR for his interest and support of this work. Thanks are also due to M. Malshe, C. Nagasubramanian, R. Narulkar, V. Karuppiah, D. Joshi, J. Ma, H. Viswanathan as well as to the Molecular Dynamics and FEM research group under Dr. Komanduri, for their support and helpful conversations. iv I wish to extend my gratitude to R. Raghav and Y. Liu for their encouragement and friendship. Special thanks are due to Dr. Gary Young for his support during this graduate study. Finally, I am indebted to my father Prof. Pandurang B. Daphalapurkar, and my mother Rupali P. Daphalapurkar for their firm belief in education and beloved concern all the time. Their encouragement and support all the time have helped me successfully fulfill this goal. v ABSTRACT A new multiscale simulation approach is introduced that couples atomistic scale simulations using molecular dynamics (MD) with continuum scale simulations using the recently developed material point method (MPM). In MPM, material continuum is represented by a finite collection of material points carrying all relevant physical characteristics, such as mass, acceleration, velocity, strain, and stress. The use of material points at the continuum level provides a natural connection with the atoms in the lattice at the atomistic scale. A hierarchical mesh refinement technique in MPM is presented to scale down the continuum level to the atomistic level, so that material points at the fine level in MPM are allowed to directly couple with the atoms in the MD. A onetoone correspondence of MD atoms and MPM points is used in the transition region, and nonlocal elastic theory is used to assure compatibility between MD and MPM regions, so that seamless coupling between MD and MPM can be accomplished. A single crystal silicon workmaterial under uniaxial tension is used in demonstrating the viability of the technique. A Tersofftype, threebody potential was used in the MD simulations. The coupled MD/MPM simulations show that silicon under nanometric tension experiences with increasing elongation in elasticity, dislocation generation and plasticity by slip, void formation and vi propagation, formation of amorphous structure, necking, and final rupture. Results are presented in terms of stress  strain relationships at several strain rates, as well as the rate dependence of uniaxial material properties. Isotropic elastoplastic constitutive model is integrated with the three dimensional formulation of MPM to simulate the material behavior of nanocrystalline materials at continuum scale. Validation of the constitutive model is done using two examples, impact load on a tensile bar and Taylor impact problem. This new multiscale computational method has potential for use in cases where a detailed atomisticlevel analysis is necessary in localized spatially separated regions whereas continuum mechanics is adequate in rest of the material. vii TABLE OF CONTENTS Chapter Page 1. INTRODUCTION………………………………………………………………… 1 1.1. Need for multiscale modeling…………………………………………….. 1 1.2. Aspects of multiscale simulation…………………………………………. 2 1.3. Overview……………………………………………………………………. 5 2. LITERATURE REVIEW…………….…………………………………………… 7 2.1. Introduction………………………………………………………………… 7 2.2. Coupled atomistic and continuum simulation approach………………. 11 2.3. Multiscale modeling of plasticity using coupled mesoplasticity and continuum simulation……………………………………………………... 38 2.4. Coupled Atomistic and Discrete Dislocation plasticity……………… .. 50 2.5. Conclusion……………………..…………………………………………. 55 3. PROBLEM STATEMENT…………….………………………………………… 57 4. COMPUTATIONAL METHODOLOGY………………………………………... 59 4.1. Introduction………………………………………………………………… 59 4.2. Molecular Dynamics (MD)………………………………………………… 60 4.3. The Material Point Method (MPM).……………….…….............………. 64 4.4. Seamless Coupling between MD/MPM.............................……………. 72 5. NUMERICAL IMPLEMENTATION OF SEAMLESS COUPLING BETWEEN MPM AND MD…………………….……………………………………………. 75 viii 6. MULTISCALE SIMULATIONS OF TENSILE TESTING………………….… 81 6.1. Introduction………………………………………………………………… 81 6.2. Coupled MD/MPM tensile testing model………………………………. 82 6.3. Deformation and failure mechanism under tensile loading…………… 84 6.4. Effects of strain rate on material properties……………………………. 92 7. ELASTOPLASTIC ANALYSIS USING MATERIAL POINT METHOD……. 96 7.1. Metal Plasticity models…………………………………………………… 96 7.2. Numerical results………………………………………………………….. 103 8. CONCLUSIONS AND FUTURE WORK ………………………………………116 BIBLIOGRAPHY......……………………………………………………………..118 ix LIST OF TABLES Table Page 1. Material parameters in Tersoff’s potential for silicon (Tersoff 1988,1989)… 62 2. Material properties of single crystal silicon (Runyan 1999)……………… 73 x LIST OF FIGURES Figure Page 1 (a). Schematic of simulations at various levels (Komanduri, 2002).……… 3 (b). Schematic of various length and time scales in tribology (Hutchings, 2002)…………………………………………………………. 3 2. Representative illustration of FE/MD handshake at the interface (Rudd and Broughton, 2000)…………………………………………………. 13 3. The stress waves propagating through the slab using a finely tuned potential energy color scale (Abraham et al. 2000)………………………………………. 14 4. Transition region for QC method (Curtin and Miller 2003)….……………. 17 5. (a) FE mesh used to model dislocationGB interaction using QC method (Shenoy et al. 1999)………………………………………………………………. 20 (b) Snapshots of atomic positions at different stages in the deformation history. (Shenoy et al. 1999)…………………………………………………. 21 6. Quasicontinuum simulation of fracture at atomic scale (Miller et al. 1998) 23 7. Nanoindentation simulation on thin film Al using QC method (Tadmor et al. 1999)……………………………………………….………… 25 8. Schematic diagram of the construction of CWM (a) from the LJ system of 1x1 units of length discretized into 256x256 points (b) from the Potts system of 2x2 units of length discretized into 128x128 points. (Frantziskonis and Deymier, 2000)…………………………………………. 31 xi 9. Connection between continuum region and atomistic region (Tan 2001).. 36 10. Nanoindentation simulation stress contour using multiscale model of plasticity. (Zbib and Rubia 2002)…………………………………………………. 43 11. Results of nanoindentation simulations using multiscale plasticity model. (a) Developed dislocation structure underneath the indenter. (b) Predicted indentation depth showing the development of surface distortions. (Zbib and Rubia 2002)……………………………………………………… 44 12. (a) CADD atomistic/continuum interface…………………………………… 52 (b) One detection element, indicating the 3 slip planes on which it lies (Curtin and Miller 2003)…………………………………………………… 52 13. (a) Nanoindentation simulation using CADD……………………………… 54 (b) Predicted forcedisplacement curve……………………………………. 54 (c) Positions of nucleated dislocations (Curtin and Miller 2002)………… 54 14. MPM grid cells with material points……………………………………….. 64 15. Hierarchical mesh refinement in 2D MPM……………………………….. 68 16. 3D MPM regular and transition cells……………………………………… 71 17. Coupled MD and MPM computational flowchart………………………… 76 18. Schematic of tensile model for coupled MD/MPM simulations………… 77 19 (a) MD atoms and MPM points initially overlapped in the transition region……………………………………………………………………… 78 (b). MPM points moving under velocity……………………………………….. 78 (c). MD atoms moving same positions as MPM points…………………….. 78 xii Figure Page 20. Animation of coupled MD/MPM model showing initial stage of tensile testing…………………………………………………………………………. 83 21. Snapshots of MD/MPM coupling simulation at strain rate 44.88 ns1…. 85 22. Snapshots of MD/MPM coupling simulation at strain rate 44.88 ns1 (magnified view)……………………………………………………………… 86 23. Snapshots of MD/MPM coupling simulation at strain rate 44.88 ns1….. 88 24. Snapshots of MD/MPM coupling simulation at strain rate 67.33 ns1….. 89 25. Snapshots of MD/MPM coupling simulation at strain rate 80.79 ns1….. 90 26. Snapshots of MD/MPM coupling simulation at strain rate 89.76 ns1…. 91 27. Engineering stressstrain curves from MD/MPM coupling simulations for silicon………………………………………………………………………… 93 28. Effects of strain rate on material properties……………………………… 94 (a). Tensile strength vs. strain rate. (b). Ultimate strain vs. strain rate. 29. 3D bar represented by collection of material points …………………….. 103 30. Applied load history…………………………………………………………. 104 31. von Mises vs time plot for elastic material model……………………….. 105 32. von Mises vs time plot for elasticperfectplastic material model with yield stress of 1.5 Pa…………………………………………………………….. 105 33. von Mises vs time plot for elasticplastic material model with power law hardening …………………………………………………………………… 105 34. von Mises stress contour at time = 10.89 s…………………………… 107 35. Stress XX contour at time = 10.89 s…………………………………… 108 xiii Figure Page 36. Displacement along X contour at time = 10.89 s……………………….. 108 37. Equivalent plastic strain contour at time = 10.89 s…………………….. 108 38. Original shape represented by material points at 0 s………………….. 110 39. Deformed shape represented by material points at 24 s……………… 110 40. Deformed shape at 76 s………………………………………………….. 110 (a). MPM points representing the deformed bar. (b). Deformed FEM mesh representing the deformed bar. 41. Equivalent plastic strain contour at 24 s…………………………………..113 (a) Axisymmetric MPM (b) Axisymmetric FEM 42. 3D Equivalent plastic strain distribution contour at material points………114 (a) At time = 24 s (b) At time = 90 s 43. Time history of the total kinetic energy…………………………………….. 115 44. Time history of the free end displacement …………………………… ….. 115 1 CHAPTER 1 INTRODUCTION 1.1. Need for multiscale simulation Material failure at all length scales experiences such processes as elastic deformation, dislocation generation and their propagation, cleavage, crack initiation and growth, and final rupture. Despite significant developments in materials simulation techniques, the goal of reliably predicting the properties of new materials in advance of material characterization and fabrication has not yet been achieved. It is also not possible to predict the properties of the material at nanoscale knowing the properties of the material at macrolevel or vice versa. This situation exists for several reasons that include a lack of reliable potentials to describe the behavior of the material at the atomistic level, computational limitations, inadequate modeling of the process, absence of scaling laws, and difficulties associated with the experimental measurement of properties even at the microscale, let alone at the nanoscale. Multiscale simulation technique can be applied to reliably predict the properties of new materials in advance of material characterization and fabrication. Scaling laws governing the mechanical behavior of materials from atomistic (nano), via mesoplastic (micro), to continuum (macro) are very 2 important to numerous DoD applications, such as the development of a new class of aircraft engine materials, or new steels for naval battle ships, or new tank armor materials for the army, or numerous microelectromechanical components for a myriad of applications, for the following three reasons: (1) information on the mechanical behavior of materials at nanolevel is not presently available as input to nanotechnology for the manufacturing of nanocomponents or microelectromechanical systems (MEMS). For example, nanostructures may possess unique properties in view of their very high surface to volume ratio, or the nanostructures might be relatively free of defects with strength approaching the theoretical values, (2) applications where two length scales of different orders of magnitude are involved. For example, one is atomistic (nano) and the other mesoplastic (micro) as in nanoindentation, or, one is mesoplastic (micro) and the other continuum (macro) as in conventional indentation, and (3) it may be possible to extend the knowledge accumulated over time on material behavior at the macro (or continuum) level to the atomistic (or nano) level, via mesoplastic (micro) level. 1.2. Aspects of multiscale simulation In recent years, combined atomistic and continuum simulations have received increasing attention due to their potential linkage between structureproperty relationships from nano to macrolevels (King el al. 1995, Phillips 1998, Shenoy et al. 1999, Smirnova et al. 1999, Smith et al. 2001, Ogata et al. 2001, and Shilkrot et al. 2002). Various computational methods were developed to 3 Figure 1 (a). Schematic of simulations at various levels (Komanduri, 2002). Figure 1 (b). Schematic of various length and time scales in tribology (Hutchings, 2002). 4 couple spatial and temporal scales. Fig. 1 (a) is a schematic of the various computer simulation methods used in the development of scaling laws as a function of spatial and temporal variables (Komanduri, 2002). At the subatomic level, we use quantum mechanics or ab initio calculations to develop the potentialenergy functions. With increase in length (or time) scale, we model material behavior using molecular dynamics (MD) and Monte Carlo (MC) simulations, then to micro or mesoplastic, and finally to continuum mechanics, such as finite element method (FEM). Fig. 1 (b) is a schematic showing various spatial and temporal scales for a typical tribological phenomenon for a wide range of applications (Hutchings, 2002). It may be noted at the outset that any scaling law developed would involve both temporal and spatial scales. At the atomistic level, MD is the preferred mode while at the continuum level, FEM is the dominant numerical technique used. Consequently, most published results on combined atomistic and continuum simulations have focused on MD/FEM coupling. At Oklahoma State University considerable work is being conducted in multiscale modeling and simulation for materials processing. The main objective of this ongoing research is to develop a fundamental understanding of indentation and tension tests at various scales and develop a computer simulation code that would bridge the gap from atomistic to continuum, via. mesoplastic or microscopic behavior. It is hoped that the outcome of this work would provide the required inputs to Computer Aided Design (CAD) software so that relevant output parameters can be used by the engineer in the design of 5 components at any level of detail, namely, large macro components, or small microcomponents, or even extremely small nanocomponents. It may be noted that the work described, herein, is part of a larger ongoing research to aid materials modeling and simulation. 1.3. Overview In Chapter 2, major multiscale modeling approaches are reviewed. The major focus of the review is on atomistic/continuum coupling. Material point method (MPM) is introduced and its potential to handle a wide range of problems is reviewed and discussed. In Chapter 4, the theoretical background and computational methodologies used in MD and MPM are reviewed. Based on these two simulation techniques, seamless coupling between these two for ability to perform coupled simulations is recognized. Chapter 5 discusses numerical implementation methodology of seamless coupling between MD/MPM. In Chapter 6, simulation approach, material properties, and results of multiscale simulations of tensile testing are presented. In tensile simulations, the effects of strain rate on the nature of deformation and material properties are investigated based on the simulated stress  strain relations for single crystal silicon. 6 Chapter 7 discusses elastoplasticity, and its computational methodology for numerical implementation in MPM. Validation results of the same are presented. Finally in Chapter 8, conclusions are made out of present investigation and offer some suggestions for future work. 7 CHAPTER 2 LITERATURE REVIEW 2.1. Introduction Continuum approach of modeling the behavior of materials, dominated the research method for numerous decades. The analysis using continuum mechanics approach implicitly assumes that, all the atomic scale dynamics and defect evolutions that determine the mechanical properties are averaged over time and large length scales. In fact, this notion stemmed from the desire to reduce the enormous number of degrees of freedom (DOF) in a given system to as minimum as possible. Thus, most of the relevant information in constitutive equations represents some averaged behavior of many, many atoms. However, the processes such as, dynamic fracture instability (Abraham et al. 1994), intersonic crack motion (Abraham and Gao, 2000), dynamic competition of dislocation emission and cleavage (Tan and Yang, 1994, Zhang and Wang, 1995), atomic inertial effect and chaotic atom motion at crack tip (Tan and Yang 1996), dislocation patterns, bifurcation phenomenon, etc., which are seen in material systems at the nano and micro scales, cannot be described by continuum theory alone. In the abovementioned processes, the length and time scales reach those of typical defects: point defects, dislocations, interfaces, and 8 grain boundaries. Hence, in case of all the abovementioned phenomena that represent material catastrophes, the averaging techniques will not yield the correct information. Nonetheless, the dawn of large scale computing is driving the art and science of modeling material phenomena into a tantalizing new direction. Instead of attempting to reduce the complexity of the material system’s behavior by a process of reduction of its DOF, one is trying to represent large DOFs, and solve for them numerically. The result is a real numerical experiment, whose outcome is not known a priori. Whether that makes sense can only be tested by confronting the outcomes of computer simulations with a limited range of experimental observations. In atomistic simulations under the given initial positions and momentum of the system, integration of Newtonian second law equation yields the total trajectory of the system. With the knowledge of the trajectories of all the atoms, one can calculate spatial and temporal distributions of energy, temperature, and pressure, as well as map out all atomistic scale processes. Advances have been made in atomistic simulation techniques, with better and more rigorous ways to approximate the quantum mechanical behavior of atoms and molecules. This improvement in connection between abinitio and MD calculations results has led to more accurate classical molecular dynamics (MD) simulations. 9 However, even atomistic analysis approach has its own length and time scale limitations. Due to explicit atomistic modeling methods, time scales are on the order of 1 fs and length scales on the order of 1 Å. Because of these very small time and length scales, typical atomistic simulations are limited to very small systems over very short times. Furthermore, there exist several multiscale methods to combine atomistic models and continuum models in a single simulation framework. These atomistic and multiscale methods have been successfully applied to investigate diverse defect structures within static or quasistatic descriptions. Moreover, even though one can model the material system with full atomic scale details, there are still some problems faced at mesoscale. As the number of atoms in a system increases, the possible local minimum energy configurations also grow rapidly (for N atom cluster the number of local minimum energy configurations grows faster than eN). Without knowing all the relative energy values of these local configurations, it is very difficult to prepare initial atomic configurations, which are most relevant to real physical systems. Secondly, the inaccuracy of an interatomic potential can introduce errors in relative energies of diverse configurations. Both of these problems seriously influence the reliability of atomistic simulations. Mesoplasticity is therefore needed to bridge this gap in between atomistic mechanisms of deformation and fracture to continuum behavior (refer Hartley 2003). Substantial progress has been made in the mesoscopic simulation area. The number of 3D dislocation loops required to represent fullscale plasticity of 10 submicron crystal is now manageable, and the solution requires integration of a few thousand equations of motion. However, the more challenging problem of polycrystalline plasticity will require additional breakthroughs, because of the conceptual difficulties of connecting dislocation dynamics (DD) to crystal plasticity models in a selfconsistent fashion. The main task is to retain the characteristic length scale from the discrete dislocation to the continuum description. Varieties of strain gradient theories have been proposed during the last decade. Most of them are formulated in a phenomenological manner and so the characteristic length is invoked into the theory, without rigorous treatment of its dislocation origins. In addition, it is also possible that the approach using a spectrum of length scales will be more appropriate to clearly establish the links between mesoscopic simulations and continuum mechanics. Recently, Ghoniem et al. (2002, 2003), Liu et al. (2004), Park and Liu (2004), and Lu and Kaxiras (2004) examined the efforts in multiscale simulation, that treat the problems involving either multiple scales in spatial or temporal or both simultaneously. Following sections give a brief review on the development of multiscale modeling and simulation of material behavior. In this review, we will be concerned with bridging length scales; while bridging time scales from the picoseconds times of atomistic vibrations to the micro, milli and larger time’s scales of various defect motions is an equally difficult but fundamentally different 11 matter. Moreover, this review mainly focuses on atomic atomistic/continuum coupling. Finally, some aspects along with recent work in coupled atomistic/continuum simulation through mesoplasticity will be reviewed in sections 2.3 and 2.4. 2.2 Coupled atomistic and continuum simulation approach Recently an alternative approach was explored that couples the atomistic and continuum approaches. The main idea behind this approach is to use atomistic modeling where the displacement field varies on an atomic scale and to use the continuum approach elsewhere. The critical part in these coupled simulations is the communication between the two descriptions of the materials. At the atomistic level, MD is the preferred mode while at the continuum level, FEM has been the dominant numerical technique used. Seamless coupling is required for the combined MD/FEM simulations, which is a key issue in multiscale simulations. The main difficulty lies in proper matching between the atomistic and the continuum regions. Since the details of the lattice vibrations, the phonons, which are an intrinsic part of the atomistic model, cannot be represented at the continuum level, conditions must be such that the phonons are not reflected at the atomisticcontinuum interface. Since the atomistic region is expected to be a very small part of the computational domain, violation of this condition quickly leads to local heating of the atomistic region and destroys the simulation. In addition, the matching at the atomisticcontinuum interface must be 12 such that largescale information is accurately transmitted in both directions. Some major multiscale atomistic/continuum simulations techniques developed are outlined below. 2.2.1 Finite element (FE) and molecular dynamics (MD) coupled simulation An early suggestion for coupling the continuum and the atomistic regimes was due to Kohlhoff et al. (1991), Hoover et al. (1992) and recently to RafiiTabar et al. (1998) and Abraham et al. (1998). In atomistic modeling of defects with long range stress fields, the boundary conditions applied to the atomistic region should be far from the defects. This results into unnecessary computations in the atomistic region far away from the defects. The other possible approach recognized was to surround the atomistic region with continuum region, which will provide flexible boundary conditions to the atomistic region. Seamless coupling is required for the combined MD/FE simulations, which is a key issue in multiscale simulations. Two methods used to combine the MD and FE region exist. One is through the imaginary surface, and the other is through the overlapping belt. The hybrid MD/FE method introduces a MD/FE imaginary surface as the handshake region to interface between the MD and the FE regions. MD atoms and FEnodes are coupled by mutual displacement boundary conditions. The kinetic and strain (potential) energy are defined for the entire system, including the handshaking interactions at the interface. The FE mesh in the MD/FE 13 handshaking region is scaled down to coincide with the atomic lattice. A onetoone correspondence of MD atoms and FEnodes is required in the interface region as shown in Fig. 2. Figure 2. Representative illustration of FE/MD handshake at the interface (Rudd and Broughton, 2000). In Fig. 2, The FE/MD interface is represented by a dashed vertical line. FE elements contributing to the overall Hamiltonian with full weight have dark shading, while those contributing with half weight have light shading. Two and threebody atomic interactions crossing the interface also carry half weight, and are shown with dotted lines. The FE mesh is gradually made coarser in the area FE region MD region MD/FE region MD/FE interface 14 far away from the interested atomistic region. Additionally, the width of the transition region is taken to be equal to the cutoff distance of the interaction potential used in the MD region. This provides a complete set of neighbors within the interaction range for all atoms in the MD region. Atoms/FEnodes that belong to the interface region not only interact via the interaction potential with the MD region but also are part of the nodes in the FE region. The positions and velocities of these atoms/FEnodes in the interface region (both on the MD side and on the FE side) must be consistent with each other. Equilibrium is guaranteed if the elastic constants in the continuum region match those of the atomistically modeled region. Figure 3 The stress waves propagating through the slab using a finely tuned potential energy color scale (Abraham et al. 2000) Kohlhoff et al. (1991) and Gumbsch et al. [1995 (a)] used the above approach in their ‘FEAt’ model for atomic scale fracture processes. They MD FE FE MD 15 obtained reasonably accurate calculation of critical energy release rates for fracture and crack tip blunting, allowing for a systematic comparison between the atomistic simulations and continuumbased criteria predicting brittle vs. ductile response. Fig. 3 shows results of coupled FE/MD simulation on rapid brittle fracture of a silicon slab flawed by a microcrack at its center and under uniaxial tension (Abraham et al. 2000). The stress wave propagation is seen in the slab through FE/MD interface with no visible reflection. Lidorikis et al. (2001) demonstrated a detailed comparison between hybrid simulation results and full MD multimillionatom simulations in the study on stress distributions in silicon/siliconnitride nanopixels. The coupled results showed a good agreement with full multimillionatom MD simulations. The atomic structures at the latticemismatched interface between amorphous silicon nitride and silicon induced inhomogeneous stress patterns in the substrate that cannot be reproduced by a continuum approach alone. However, the approach by Kohlhoff et al. (1991) does not support continuum dislocations in the finite element region, nor can it easily adapt to the size of the atomistic region during the simulation. The FEAt model requires that the atomistic region be specified at the start of the computation, thereby limiting any nonlinear deformation (including dislocations) to a predefined region The other method, which uses an overlapping belt, was found to be more flexible. In the overlapping area each node is surrounded by a group of atoms. 16 The forces exerted on the atoms in that region, due to the interaction with the MD region, make up the external forces on the FEnodes for the FE simulation. Tan and Yang (1994) and Yang et al. (1994) developed a combined continuum and atomistic simulation approach to simulate the transmission of crack tip dislocations from the atomistic assembly to the overlapping continuum. Noguchi and Furuya (1997) proposed a method that combines MD at the crack tip with a linear elastic outer domain in simulating elasticplastic crack advance. MD region was successfully moved as the crack propagated (Furuya and Noguchi, 1998, 2001). 2.2.2 Quasicontinuum method. Quasicontinuum (QC) theory developed by Tadmor et al. (1996a, 1996b) has been applied to a wide range of problems including dislocation structures, interactions of cracks with grain boundaries (Shenoy et al. 1998), nanoindentations (Shenoy et al. 2000, Smith et al. 2000, Smith et al. 2001), dislocation junctions (Rodney and Phillips 1999), atomistic scale fracture process (Miller et al. 1998). The results obtained are quite encouraging. Ortiz and Phillips (1999) reviewed the QC theory. Ortiz et al. (2001) established the same for multiscale simulations. QC seamlessly couples atomistic and continuum regions. In this approach the atomistic and finite element calculations are performed concurrently rather than in sequence, and therefore it is categorized as a concurrent multiscale approach. The QC theory initializes from a fundamental conventional atomistic model and systematically eliminates redundant degrees of 17 freedom by introducing kinematic constrains. The simulation initializes from atomistic scale and coarsens to describe the continuum region. While doing this, kinematic constraints are thoughtfully introduced by preserving the full atomistic description wherever required. Figure 4. Transition region for QC method (Curtin and Miller 2003) The QC method defines two types of atoms, ‘local representative atoms’ and ‘nonlocal representative atoms’, rather than identifying atomistic and continuum regions. In practice, however, the regions containing nonlocal representative atoms are essentially equivalent to the fully atomistic regions of other methods. Similarly, a local representative atom is coincident with either a continuum FE node or an atomic position near a Gauss point used to define the 18 energy of the continuum element. The language of nonlocal and local clearly associates the nonlocal atoms with ‘real’ atoms and the local atoms with the local FE region. The transition region between nonlocal and local representative atoms in the QC model is depicted in Fig. 4. The interface atoms are included in the atomistic energy but are also FE nodes of the continuum. The pad atom positions are dictated by interpolation from the FE nodal positions. The pad atom energies are not included but the energies of the real and interface atoms include their interactions with the pad atoms. To avoid overcounting of the energy of the interface atoms/nodes, the energies of the continuum elements adjacent to the interface (i.e. elements that have nodes corresponding to the interface atoms as shown in grey in Fig. 4) are weighted differently in the total potential energy sum. The total potential energy of the QC model is then obtained by summing the energies of all atoms in the atomistic region and at the interface and all elements in the continuum domain In this method, interatomic interactions are used in the model via the discrete calculation based on a local state of deformation. It starts by computing the energy of solid as a function of atomic positions. The configurations space of the solid is reduced to a subset of representative atoms. The over all equilibrium equations are then obtained by minimizing the potential energy described for the reduced configuration space. The total energy, in principle, needs to be calculated by summing the energy over all the atoms in system. N i 1 tot i E E , (2.1) 19 where N is the total number of atoms in the body. This full summation is bypassed by introducing approximate summation rules. The lattice quadrature analog of Eqn. (2.1) is written as Nr i 1 tot i i E n E where, ni is the quadrature weight that signifies how many atoms a given representative atom stands for in the description of the total energy, and Ei is the energy of ith representative atom which is summed over Nr representative atoms to get the total energy. In QC approach, displacement fields are obtained using FE method. The positions of the representative atoms are also calculated. From the new displacement field, again, energies are determined using atomistic calculations and the simulation continues. The method has been applied to single crystal FCC metals and shown to reproduce lattice statics results for a variety of line and surface defects (Tadmor et al. 1996b). Shenoy at al. (1998,1999) generalized the QC method and extended it to treat polycrystalline materials. As the first test of the model, they computed the grain boundary (GB) energy and atomic structures for various symmetric tilt GB’s in Au, Al, and Cu using both direct atomistic calculations and the model calculations. They found excellent agreement between the two sets of calculations, indicating the reliability of the model for their purpose. The generalized method was applied to a number of problems including the interaction of dislocations with grain boundaries in Al and the mechanics of steps on grain boundaries (Shenoy et al. 1998). In the study of Al, they used 20 nanoindentation induced dislocations to probe the interaction between dislocations and GB’s. Figure 5(a) FE mesh used to model dislocationGB interaction using QC method (Shenoy et al. 1999) Specifically, they considered a block oriented so that the (111) planes are positioned to allow for the emergence of dislocations which then travel to the GB located at ~200Å beneath the surface. Fig 5 (a) shows FE mesh for nanoindentation model used for multiscale simulation using QC method. The surface marked AB is rigidly indented in order to generate dislocations at A. First, the energy minimization is performed to obtain the equilibrium configuration of the GB, and then the FE mesh is constructed accordingly. Full atomistic description is used for the region that is expected to participate in the dislocation 21 GB interaction, while in the far fields the mesh is coarser. The displacement boundary conditions at the indentation surface are then imposed onto this model structure, and after the critical displacement level is reached, dislocations are nucleated at the surface. They use the EAM potential to calculate the energies at atomistic level. They observed closely spaced (15 Å) Shockley partials nucleated at the free surface. Figure 5(b) Snapshots of atomic positions at different stages in the deformation history. (Shenoy et al. 1999) As seen in [Fig 5(b)], the partials are subsequently absorbed at the GB with the creation of a step at the GB and no evidence of slip transmission into the adjacent grain is observed. The resultant structure can be understood based on the concept of displacement shift complete lattice associated with this symmetric tilt GB. As the load is increased, the second pair of Shockley partials is 22 nucleated. These partials are not immediately absorbed into the GB, but instead form a pileup [Fig. 5(b)]. The dislocations are not absorbed until a much higher load level is reached. Even after the second set of partial dislocations is absorbed at the GB, there is no evidence of slip transmission into the adjacent grain, although the GB becomes much less ordered. The authors argued that their results give hints for the general mechanism governing the absorption and transmission of dislocations at GB’s. In the same work they also studied the interaction between a brittle crack and a GB and observed stressinduced GB motion (the combination of GB sliding and migration). In this example, the reduction in the computational effort associated with the quasicontinuum thinning of degrees of freedom is significant. For example, the number of degrees of freedom associated with the mesh is about 104, which is three orders of magnitude smaller than what would be required by a full atomistic simulation (107 degrees of freedom). Miller et al. (1998) extended the method to study fracture mechanics at the atomic scale. Fig. 6 shows a crack tip subjected to mode I loading approaching a grain boundary in a nickel bicrystal, with increasing levels of external load. The colors indicate levels of slip from zero (red) to maximum (blue). The snapshots show initial dislocation from the grain boundary (t=3), followed by crack extension by cleavage (t=4) and in the end grain boundary migration to the crack tip (t=6). They also investigated the effect of grain orientation on fracture and the interaction of cracks with grain boundaries. 23 Rodney and Phillips (1999) extended the method to three dimensions and used it to study junction formation and destruction between interacting dislocations. Further, the method was also extended to treat complex Bravais lattice materials (Tadmor et al. 1999). Figure 6 Quasicontinuum simulation of fracture at atomic scale (Miller et al. 1998) The authors refer to using the “local” regime of the QC method, which refers to the case where each FE node represents a very large number of 24 atomistic degrees of freedom, which is modeled as corresponding to an infinite solid homogeneously deformed according to the local strain at the node. In this limit, it is advantageous to employ effective Hamiltonians to compute energetics for the representative atoms. Such Hamiltonians can be constructed from abinitio calculations, and may include physics that atomistic simulations based on classical interatomic potentials (such as EAM) are not able to capture. For example, by constructing, an effective Hamiltonian parameterized from ab initio calculations, Tadmor and his group was able to study polarization switching in piezoelectric material PbTiO3 in the context of the quasicontinuum model in the local limit (Tadmor et al. 2002). With this effective Hamiltonian, it was shown that the behavior of a largestrain ferroelectric actuator under the application of competing external stress and electric fields could be modeled successfully, reproducing, for example, all the important features of the experimental strain vs. electric field curve for the actuator. The advantage of these simulations is that they also provide insight into the microscopic mechanisms responsible for the macroscopic behavior, making possible the improvement of design of the technologically useful materials. QC was also used to study nanoindentation into silicon single crystals (Smith et al. 2000, 2001) and of aluminum thin film (Tadmor et al. 1999). Fig. 7 shows nanoindentation simulation of a thin film aluminum that was carried out. Dislocations are seen emitted from indenter corners (indenter not in the picture). 25 Figure 7 Nanoindentation simulation on thin film Al using QC method (Tadmor et al. 1999) Knap and Ortiz (2001) presented a streamlined and fully nonlocal threedimensional version of the QCtheory. However in spit of all this, the QC method does not allow for the continuum description of a dislocation and therefore every dislocation in the QC model requires a full complement of atomistic degrees of freedom around its core and slip plane. This means that only a few dislocations can be present in a QC simulation before the computational cost becomes comparable to that of the fully atomistic description. Additionally the computational effort in the QC approach quickly approaches that of the fully atomistic model once a relatively small number (less than 100) of dislocations are developed in the microstructure. 26 2.2.3 Lattice Green’s function. Thomson et al. (1992), Zhou et al. (1993) and SchiØtz et al. (1997) used this multiscale approach for calculating the lattice defects. In this approach, the actual number of atomic degrees of freedom is not changed. They treat a large portion of atomic response as linear instead of nonlinear while they are still explicitly modeled. As a result they considerably reduced the computational time required. Rao et al. (1999) extended the Green’s function technique to two and three dimensions. They demonstrated this technique by applying it to relax the boundary forces in the simulations of crossslipped core structures of (a/2)[110] screw dislocations of FCC structure. Further, MasudaJindo et al. (2001) used this approach to study the mechanical properties of nanoscale materials. The initial atomic structures of the dislocations and cracks were determined both from the elastic solutions as well as from those by lattice Green’s function method for the infinite systems. They calculated the Green function for the defective lattice, with dislocation and crack, by solving the Dyson equation, appropriate for absolute zero temperature. The thermal expansion and the temperature dependence of the interatomic force constants are determined using the statistical moment method and they are taken into account in the lattice Green’s functions. With this methodology, they investigated the strength and fracture properties for the nanocrystalline materials like semiconductor quantum wire and carbon related materials like graphenes and nanotubes. The O(N) tightbinding molecular dynamics (TBMD) method is 27 used to analyze the reconstruction of atomic bonding near the crack tip as well as the cleaved surface. . Overall, Lattice Green’s function approach is used for static systems and is effective in seeking longterm balanced configurations. 2.2.4 Coarsegrained molecular dynamics. Rudd and Broughton (1998) took into account the limitation that the continuum based origins of the FE method lead to difficulties in a smooth transition from the atomistic to the continuum. Thus, they reformulated the continuum region using a method called coarsegrain molecular dynamics (CGMD), in which they developed the degreeoffreedom reduction as a more natural extension of the underlying discrete molecular dynamics. CGMD substitutes FE mesh used outside the atomistic region, and it connects to MD at the atomistic scale. In this approach, the equations of motion do not rely on the assumptions of continuum elasticity theory. As a result, CGMD was found more appropriate for mesoscopic elastic solids. The model aims at constructing scaledependent constitutive equations for different regions in a material. In their methodology, the material of interest is partitioned into cells, whose size varies so that in important regions a mesh node is assigned to each equilibrium atomic position, whereas in other regions the cells contain many atoms and the nodes need not coincide with atomic sites. The CGMD approach produces equations of motion for a mean displacement field defined at the nodes by defining a conserved energy functional for the coarsegrained system as a 28 constrained ensemble average of the atomistic energy under fixed thermodynamic conditions. The key point of this effective model is that the equations of motion for the nodal (mean) fields are not derived from the continuum model, but from the underlying atomistic model. The nodal fields represent the average properties of the corresponding atoms, and equations of motion (in this particular case Hamilton’s equations) are constructed to describe the mean behavior of underlying atoms that have been integrated out. One important principle of CGMD is that the classical ensemble must obey the constraint that the position and momenta of atoms are consistent with the mean displacement and momentum fields. When the mesh nodes and the atomic sites are identical, the CGMD equations of motion agree with the atomistic equations of motion. As the mesh size increases some shortwavelength degrees of freedom are not supported by the coarse mesh. However, these degrees of freedom are not neglected entirely, because their thermodynamic average effect has been retained. This approximation is expected to be good if the system is initially in thermal equilibrium, and the missing degrees of freedom only produce adiabatic changes to the system. The Hamiltonian was derived originally for a monoatomic harmonic solid, but can be easily generalized to polyatomic solids (Rudd and Broughton 1998). After deriving the equations of motion from the assumed Hamiltonian for a particular system, one can perform the CGMD for the nodal points. As a proof of principle, the method was applied to onedimensional chains of atoms with periodic boundary conditions where it was shown that the CGMD 29 gives better results for the phonon spectrum of the model system compared to two different FE methods (Rudd and Broughton 1998). A variety of other calculations have also been performed with the CGMD to validate its effectiveness (Rudd and Broughton 2000, unpublished). Rudd (2001) developed a nonconservative version of CGMD to remove unphysical elastic wave reflection. Further, Rudd (2002) used CGMD to study the behavior of submicron MicroElectroMechanical Systems (MEMS), especially microresonators, often called NanoElectroMechanical Systems (NEMS). Although the CGMD has proven to be reliable in the description of lattice statics and dynamics, the implementation of the model varies from system to system. This is because different approximations have to be made to the Hamiltonian that represents a particular system. On the other hand, such approximations can be estimated and controlled in the CGMD method. This advantage makes the CGMD method a good substitute for the FE method in the MAAD approach when a high quality of FE/MD coupling is required. The CGMD approach resembles the quasicontinuum model in the sense that both approaches adopt an effective field model, an idea rooted in the renormalization group theory for critical phenomena. In both approaches, less important (long wavelength) degrees of freedom are removed while the effective Hamiltonian is derived from the underlying finescale (atomistic) model. Although both approaches are developed to couple FE and atomistic models, the quasicontinuum method is mainly applicable to zero temperature calculations while the CGMD is designed for finite temperature dynamics. The 30 quasicontinuum model has shown its success in many applications, but the CGMD approach has yet to show its wider applicability and versatility. The CGMD results obtained were limited to validation problems that demonstrated the capabilities of the method. 2.2.5 Compound wavelet matrix method. This method was found suitable for multiscale investigations. Frantziskonis and Deymier (2000) designed compound wavelet matrix (CWM) method to bridge multiscale models over different ranges of spatial and temporal scales. They applied CWM to analyze the microstructure of 2D polycrystalline systems as well as their evolution during normal grain growth. On the finer scale MD simulation with LennardJones (LJ) system, and on the coarser scale MC simulation of the Qstates Potts model was used. The basic idea used was to produce matrices of the wavelet coefficient from energy maps representing the spatial distribution of the local excess energy in the microstructures obtained with both methodologies. The full description of the material was then obtained by merging the matrices of the wavelet coefficients representing the material at different scales through CWM method. The CWM then characterizes the materials over a range of scales that is the union of the scales treated by the two methodologies. This method possesses several advantages. First, it does not assume a priori that a collection of small microscale systems is equivalent to a microscalebased model of a large system. 31 Second, the simulation time of the coarsest methodology is not controlled by the methodology with the slowest dynamics. Figure 8 Schematic diagram of the construction of CWM (a) from the LJ system of 1x1 units of length discretized into 256x256 points (b) from the Potts system of 2x2 units of length discretized into 128x128 points. (Frantziskonis and Deymier, 2000) To present an illustrative example of the CWM method applied to 2D grain growth, an MD simulation of a 2D LJ system and a MC simulation of a Qstates Potts model that can overlap over a range of spatial and time scales were presented, Fig. 8. The bridging of these simulations was accomplished by the 32 overlap in scale of the ‘mesoscopic’ Qstates Potts model with the atomistic LJ model and constructing the CWM. In Fig. 8 the arrows indicate the substitution (transfer of matrix statistics) of matrices in forming the compound matrix. Only three such substitutions are shown for illustration (Frantziskonis and Deymier, 2000). 2.2.6 Other approaches. Another recent development in coupled method includes algorithms to pass defects between that atomistic and continuum regions developed by Shilkrot et al. 2002 and given the name, the coupled atomistics and discrete dislocation (CADD) method. The CADD approach is similar to the FEAt approach, in that it begins from a separation of the physical problem into spatial regions that are modeled either by a fully atomistic description or by continuum finite elements. CADD extends these methods, however, by allowing for the presence and movement of continuum discrete dislocations in the continuum regime that interact with each other and with the atoms in the atomistic region via their elastic fields. The philosophy behind the CADD method is that atomistic resolution is required in some regions of a material, but away from those regions the deformation is well described by elasticity plus the plasticity due to the motion of welldefined continuum dislocations whose interactions can be described by elastic continuum fields. 33 Some other concurrent approaches are similar to the MAAD method but concentrate on linking two different length scales rather than three. For example, Bernstein and Hess (2001) have simulated brittle fracture of Si by dynamically coupling empiricalpotential MD and semi classical TB approaches. 2.2.7 Lattice material point method (LMPM) 2.2.7.1 Material Point Method (MPM). A new computational method, namely, the material point method (MPM) was developed by Sulsky and others (1995,1996) at the University of New Mexico with support from the Sandia National Laboratories. MPM was developed from FLIP (Fluid  implicit  particle) particle – in – cell (PIC) method by extending the computational fluid dynamics capability of the code, to solving solid mechanics problems. Sulsky et al. (1994), initially considered application of MPM to 2 dimensional impact problems and demonstrated the potential of MPM. Sulsky and Schreyer (1996) gave a more general description of MPM, along with special considerations relevant to axisymmetric problems. The method utilizes a material or Lagrangian mesh defined on the body under investigation, and a spatial or Eulerian mesh defined over the computational domain. The set of material points making up the material mesh is tracked throughout the deformation history of the body and these points carry with them a representation of the solution in a 34 Lagrangian frame. Interactions among these material points are computed by projecting information they carry onto a background finite element mesh where equations of motion are solved. They reported that the material point method does not exhibit locking or an overly stiff response in simulations of upsetting. In MPM, a solid body is discretized into a collection of material points, which together represent the body. As the dynamic analysis proceeds, the solution is tracked on the material points by updating all required properties, such as position, velocity, acceleration, stress, and temperature. At each time step, the material point information is extrapolated to a background grid, which serves as a computational scratch pad to solve the equations of motions. The solution on the gridnodes is transferred using finite element shape functions, on the material points, to update their values. This combination of Lagrangian and Eulerian method has proven useful for solving solid mechanics problems involving materials with history dependent properties such as plasticity or viscoelasticity effects. MPM is amendable to parallel computation (Parker 2002), implicit integration method for unconditionally stable time increments (Guilkey and Weiss 2003, Sulsky and Kaul 2003, Cummins and Brackbill 2001) and alternative interpolation schemes using C1 continuous interpolating functions for improved smoothness, without numerical noise, in field properties (Bardenhagen and Kobar 2000). The computational methodologies of material point method as well as of molecular dynamics are discussed in Chapter 4. 35 2.2.7.2 Coupled MPM/MD coupling approach Tan (2001) proposed this method. The basis for LMPM, coupling MD and MPM, is that a continuum material point is also an aggregate of atoms. MD/MPM connection was realized by the Lattice Material Point (LMP) and the background grid. LMP has two viewpoints  a lattice from an atomistic perspective and a material point from a continuum perspective. Providing a smooth transition or handshake between these two representations is the spirit of LMPM. From an atomistic point of view, a LMP represents a lattice, a set of atoms. Each LMP has corner points to define the occupied space. The corner point updates the position according to the interpolated velocity field from the MPM grid. The LMP continuum stresses are the average of atomic stresses, which are evaluated through statistical average of the force on each atom from the local neighboring atoms. LMPM uses a uniform background grid that hosts both continuum cells and atomistic cells. The simulation region is divided into two regions  an atomistic region and a continuum region; grid nodes connect the two regions. Nodes at the interface between the continuum region and atomistic region, such as node m in Fig. 9, collect stress and temperature information from LMP in both continuum and atomic sides. On the continuum side, the continuum LMP (gray points) interpolates the value from grid nodes and updates the values at material points (black points). At atomic side, the MD simulation needs the neighboring atomic position information to compute the potential and the forces on each atom (black points). 36 Figure 9. Connection between continuum region and atomistic region (Tan 2001). Using LMPM Tan (2001) simulated an intergranular crack that propagates to an intergranular SiC under mode I extensional load. LMP was used to convert atomistic dislocations to continuum ones. When an atomistic dislocation moves to the continuum region, it was substituted in situ by a dislocation of the same Burgers vector but embedded in the MPM described continuum. As the consequence of the newly created continuum dislocation, the extra halfarray of atoms in the atomistic assembly were removed from the MD calculation to retain the global conservation of mass. In the continuum region, dislocation motion was directed by the crystal slip system with speed given by the dislocation dynamics curve. 37 Ayton et al. (2001) developed a computational methodology, which included the continuum level simulation using MPM and the microscopic method of nonequilibrium MD. They examined the behavior of a membrane bilayer/solvent system using this microtomacro dynamical feedback simulation technique. With this method, two simulations at different time and length scales were interfaced into a unified simulation methodology. The interface was accomplished via an information transfer where selected material properties (transport coefficients) and state parameters (density) were calculated in one spatial/temporal regime and then used as initial input in another leading to a closed feedback loop. The multiscale methods based on Lattice Green’s function and QC methods are static methods. These give the longterm balanced configurations. Also time scale is not an issue in these approaches. However dynamic effects are important at nanoscale, (refer Tan 2003). Dynamic analysis at nanoscale is possible using coupled FE/MD and LMPM methods. The material point/atoms linkage in LMPM has some similarity with CGMD, where a material point can be looked as a coarsegrained atom group. However, in FE/MD coupling, one—toone correspondence of FEnode and MDatom provides extra constraint and blocks the moving of atomic defects into continuum region (Tan 2001). Even in CGMD this may become a difficult issue. Additionally, under large deformations, the FE elements become distorted and it is difficult to implement the entire computational process in the FE simulations. 38 2.3 Multiscale modeling of plasticity using coupled mesoplasticity and continuum simulation. Two complimentary approaches in the area of mesoscale analysis have been advanced. First of the two approaches involves dislocation dynamics and the second approach is using continuously dislocated continuum theory or statistical mechanics approach. Recently numerous efforts have been devoted for coupled mesoscopic and continuum scale simulations. Herein the work contributed by Zbib and Rubia (2002), Hartley (2001, 2003), Khan et al. (2004),, Abu AlRub and Voyiadjis (2004), LeSar and his coworkers is reviewed. Dislocation Dynamics (DD) has been proven useful for simulating material behavior at mesoscale. In this strategy, instead of simulating the dynamics of atomic systems, the dynamics of defect ensembles are simulated in the material. Examples of this approach are the dynamic simulations of interacting cracks in brittle materials, or dislocations in crystalline materials. DD theory involves incorporating nonlinear interaction among dislocations and interaction of dislocations with interfaces. Lepinoux and Kubin (1987) were the first to develop twodimensional models of dislocations, although an alternate version. More developed versions were then proposed (Ghoniem and Amodeo 1988a,b, Amodeo and Ghoniem 1990a,b, Gulluoglu et al. 1989, Lubarda et al. 1993, Barts and Carlsson 1995, Wang and LeSar 1995). During this time the modeling 39 capabilities were limited to twodimensional and consisted of infinitely long and straight models of dislocation in an isotropic infinite elastic medium. Using DD, complex phenomena such as dislocation cells, slip bands, microshear bands, persistent slip bands and dislocation tangles, were proved to dictate the material properties. To understand more realistic features of the microstructure in crystalline solids, Kubin and coworkers (Kubin and Canova 1992, DeVincre and Kubin 1994), pioneered the development of 3D lattice dislocation dynamics. More recent advances in this area have contributed to its rapid development (Zbib et al. 1998, Ghoniem et al. 2000, 2001). Another complimentary approach for mesoscale analysis has been based on a continuously dislocated continuum theory or the statistical mechanics approach (refer Needleman 2001) and is also termed as mesoplasticity. In this approach, models involving interactions of dislocation ensembles and the manner they influence material properties are included in the evolution equations based on classical continuum mechanics framework using representative volume element (RVE). These governing equations along with boundary conditions are solved for the complete description of the deformation problem. The mesoplastic constitutive relation was first proposed by Taylor (1938). Kröner (1970) proposed a theory of global plastic multiaxial deformation that incorporated a description of the local dislocation structure in the form of state variables. Mesoplastic constitutive relation has been refined and improved by many researchers, involving Hill and Rice (1972), and Asaro and Needleman (1985). With the 40 increasing computational power, the relation can be implemented in some FEM codes to solve complex threedimensional problems. For example, using the mesoplastic theory, Peirce et al. (1982) analyzed numerically the nonuniform and localized deformations of ductile single crystals subjected to tensile loading. Huang (1991) developed a usermaterial subroutine incorporating mesoplasticity in the commercially available ABAQUS FE implicit program. This twodimensional subroutine has implemented the theoretical framework of Hill and Rice (1972). Kalidindi et al. (1992) developed an implicit timeintegration procedure mainly based on the constitutive model of Asaro and Needleman (1985). Kalidindi et al. simulated the loaddisplacement curve in a planestrain forging on copper that captured successfully major features of the experimental data. They reported jumps on the calculated curve although no jumps were observed from the experimental data. Yoshino et al. (2001) applied the mesoplasticity theory in 2D FEM to simulate the dislocation generation and propagation during indentation of single crystal silicon. Kysar (2001) studied the crack growth along a copper/sapphire bicrystal interface. Wang et al. (2004) studied the dependence of nanoindentation pileup patterns and microtextures on the crystallographic orientation on single crystal copper using a conical indenter. They, however, reported an order of magnitude’s deviation of the numerical loaddisplacement results from experimental data. Fivel et al. (1998) developed a 3D model to combine discrete dislocation with FEM for nanoindentation simulation. 41 Above mentioned models are based on basic deformation mechanisms and aim at predicting strength properties of metals at small length scales. However sizedependent problems which include, modeling of dislocation boundaries, dislocation in heterogeneous structures like multiple layer and thin film structures, dislocation interaction with interfaces and associated shape changes with lattice rotations, and additionally deformation in nanostructured materials, localized deformation and shear bands still cannot be addressed using the technique of DD alone due to computational limitations. Hence coupled DD and continuum scale approach is being developed. The main difficulty in these models, among other things, lies in connecting the observed parameters (used in continuum models), such as stress and strain, to the collective behavior of dislocation groups and their interaction with particles and interfaces. While within the continuum mechanics framework, the governing equations of the material response are developed based on a representative volume element (RVE) over which the deformation field is assumed homogeneous, the DD models are presumably capable to describing the heterogeneous nature of the deformation field within the RVE. Thus with a proper homogenization theory, or field averaging, one can couple the discrete DD models with continuum approaches and provide a rigorous framework for analyzing deformations at small scales, where surfaces and interfaces are of great importance in determining material response. van der Giessen and Needleman (1995) developed this approach. They coupled a twodimensional discrete dislocation dynamic model with continuum finite element model. Further, 42 Zbib and Rubia (2002) extended the framework to merge the two scales, namely nanomicroscale, where plasticity is determined by explicit threedimensional dislocation dynamics (DD) analysis, and the continuum scale, where energy transport is based on basic continuum mechanics laws. The output of this approach is a coupled threedimensional discrete DD with continuum elastoviscoplastic model. A method based on superposition principle is adopted for modeling dislocation surface interactions in both homogeneous and heterogeneous materials. Finite element method (FEM) is the technique adopted for continuum scale analysis. Three models are combined together, the DD model, the solid mechanics and the heat transfer model. Their hybrid multiscale model of plasticity is capable of explicitly modeling complex smallscale plasticity phenomena, including surface effects, ledges, microshear bands, dislocation boundaries, deformation of thin layers. Their results on multiscale nanoindentation simulation, for investigating deformation and strength in thinlayered structures are shown in Fig. 10 and 11. The test sample used was Fe3% Si single crystal and Oxide thin film of thickness 10.0 nm on it and the load was applied using a Berkovich indenter (192 nm tip radius of curvature) and in a direction normal to (001) plane. For more details, refer to Zbib and Rubia (2002). 43 Figure 10 Nanoindentation simulation stress contour using multiscale model of plasticity. (Zbib and Rubia 2002) The stress field using multiscale nanoindentation simulation is shown Fig 10. For the case of random distribution of dislocation sources the evolution of dislocations is shown in Fig 11 (a) and the depth during nanoindentation as predicted by their hybrid multiscale model is shown in Fig 11(b). The authors comment that the resulting dislocation structure and predicted depth are consistent with the experimental observations reported by Zielinski et al. (1995) and Bahr and Zbib (unpublished). 44 Figure 11 Results of nanoindentation simulations using multiscale plasticity model. (a) developed dislocation structure underneath the indenter. (b) Predicted indentation depth showing the development of surface distortions. (Zbib and Rubia 2002). Further, Khan, Zbib and Hughes used multiscale DD plasticity to model planar dislocation boundaries in deformed FCC single crystals. For multiscale simulations, they constructed a pure tilt boundary and experimentally observed extended geometrically necessary boundaries (GNBs) within the representative volume element (RVE). Coupling between discrete DD analysis approach with 45 continuum finite element is obtained to correct for the boundary conditions and image stress. Their approach is that, they calculate the PeachKoehler force acting on a dislocation segment from the stress field of all other dislocations and the applied stresses. And this is calculated while the plastic deformation of the single crystal is obtained explicitly accounting for the evolution of a multitude of dislocation loops and curves. Then following the standard finite element procedure and using linear interpolation shape functions over the segment, the PeachKoehler force per unit length is distributed equally to the nodes. Thus, once all the forces are assembled, the net force on each node would have contributions from all the segments connected to it. The Dislocation nodes move simultaneously in the glide direction over a characteristic time corresponding to the least time increment required for an interaction to take place. Since the governing ewuation of glide motion for each dislocation node is nonlinear, the motion and interaction of an ensemble of dislocations in a 3D crystal is integrated over time yielding macroscopic plastic distortion. The authors conclude that the right boundary condition of the RVE is critical in modeling GNBs and their longrange stresses. They also considered the effects of various numerical factors such as domain length and mesh sensitivity. Additionally the effect of changing the spacing between two dislocation boundaries on the selfstress field and the stability, particularly in the space between the two dislocation boundaries, is also discussed. 46 Towards an effort to couple the classical continuum mechanics with mesoplasticity, Hartley (1985, 2001) proposed the concept of a dislocation mobility tensor. Recently Hartley (2003) developed a method for linking thermally activated dislocation mechanisms of yielding with continuously dislocated continuum plasticity theory as developed by Kröner (1958) and Bilby (1955). For information transfer in coupled simulations, he proposed that the total distortion needs to be decomposed into components that are due to dislocation motion alone and to the subsequent lattice distortion caused by the dislocations and external boundary conditions. The volumetric distortion associated with movement of dislocations should be used to connect to continuum theories. To attain the same, he introduced two quantities, the dislocation mobility tensor and the dislocation velocity vector, that characterizes the geometrical nature and stress dependence of dislocation ensembles in a way that permits their inclusions in continuum theories of elastoplastic deformation. The magnitude of the dislocation density vector for a slip system describes the net length, or volume average, of dislocation lines having the same Burgers vector in a RVE. The orientation of the vector is determined by the relative proportion of the total length due to screw and edge components. The vector is constructed by first defining a dislocation distribution vector for an arbitrary direction as the product of a unit vector in the direction and the total length per unit volume of dislocation segments parallel to that direction. Components of the dislocation density vector relative to a reference coordinate system are obtained by integrating the corresponding components of the dislocation distribution vector over all possible 47 orientations. The magnitude of the screw component and the orientation of the vector relative to the slip direction are determined by projecting the dislocation density vector parallel to the slip direction. Defining the dislocation density vector in this manner permits the definition of the net Peach–Koehler force per unit volume as the virtual force on the dislocation configuration due to internal and external sources of stress. It is important to note that the dislocation density vector defined in this fashion only describes the screw–edge character of the net dislocation population. It does not contain details of the local geometry of the dislocation arrangement, such as the relative positions of dislocation segments; hence, it cannot be used for calculations that depend on these values. Above concepts can be applied directly to the kinematics and kinetics of dislocation dynamics calculations that produce detailed pictures of the dislocation arrangement in a computational volume. The description of dislocation density can be employed to insert the results of these calculations into continuum descriptions of macroscopic material behaviour in order that the dislocation dynamics calculations can be employed to derive constitutive equations applicable to macroscopic structures. The dislocation mobility tensor provides a convenient form for the inclusion of specific dislocationbased microscopic constitutive relations into a model that can then be used as the basis for dislocation dynamics calculations or further analytical treatments. Constitutive equations that describe the relationship of the applied stress to the associated dislocation strain rate at the onset of observable permanent deformation can be developed in terms of internal structure of the material from a model based on 48 thermal activation of dislocations over barriers in their glide plane. Using this approach Hartley demonstrate how stress, temperature and material parameters enter into an expression for the flow law using the intersection of a gliding dislocation with a forest dislocation as an illustration. Recently, Abu AlRub and Voyiadjis (2004) made use of micro to nanoindentation experiments to determine analytically the material intrinsic length scale of strain gradient plasticity theory. The enhanced gradient plasticity theories formulate a constitutive framework on the continuum level that is used to bridge the gap between the micromechanical plasticity and the classical continuum plasticity. They are successful in explaining the size effects encountered in many micro and nanoadvanced technologies due to the incorporation of an intrinsic material length parameter into the constitutive modeling. However, the full utility of the gradient type theories hinges on one’s ability to determine the intrinsic material length that scales with strain gradients, and this study aims at addressing and remedying this situation. Based on the Taylor’s hardening law, a micromechanical model that assesses a nonlinear coupling between the statistically stored dislocations (SSDs) and geometrically necessary dislocations (GNDs) is used in order to derive an analytical form for the deformationgradientrelated intrinsic lengthscale parameter in terms of measurable microstructural physical parameters. They also present a method for identifying the lengthscale parameter from micro and nanoindentation experiments using both spherical and pyramidal indenters. The deviation of the 49 Nix and Gao (1998) and Swadener et al. (2002) indentation size effect (ISE) models’ predictions, from hardness results at small depths for the case of conical indenters and at small diameters for the case of spherical indenters, respectively, is largely corrected by incorporating an interaction coefficient. It thus compensates for the proper coupling between the SSDs and GNDs during indentation. They also present experimental results which show that the ISE for pyramidal and spherical indenters can be correlated successfully by using the proposed model. LeSar and coworkers (LeSar and Rickman 2003) have devoted some efforts in developing coarsegraining strategies due to sizelimitation imposed by direct DD simulations. They have begun to tackle this issue by investigations of a number of specific problems, namely the combined motion of dislocations and solutes, and the problem of the longrange interactions of blocks of dislocation ensembles. In the limit where the solute mobility is high relative to that of dislocations, dislocation–solute atmosphere pairs can be treated as quasiparticles, and the degrees of freedom associated with individual solute atoms can be eliminated (Rickman et al. 2003). Thus, an effective mobility can be derived and can be related to the details of such interaction. However, in the other extreme limit of low solute mobility compared with dislocations, the solute atoms function essentially as static traps for dislocations. In this case, LeSar and coworkers modeled slip as a discrete onedimensional ‘birthanddeath’ stochastic process. These examples show that it may be possible to average out fast time 50 variables in DD simulations. The issue of spatial averaging has been recently dealt with by LeSar and Rickman (2002) and Wang et al. (2003b). Starting from the work of Kossevich (1979) on the interaction energy of systems of dislocations, an energy expression in terms of the dislocation density tensor (Rickman et al. 2003) is derived. The basic idea in this approach is to divide space into small averaging volumes to calculate the local multipole moments of the dislocation microstructure and then to write the energy (LeSar and Rickman 2002) or the force (Wang et al. 2003b) as an expansion over the multipoles. Numerical implementation of the multipole expansion for the Peach–Kohler force on a test dislocation separated from a volume containing a dense dislocation aggregate revealed that such a technique is highly accurate and can result in significant coarse graining in current DD simulations. 2.4 Coupled Atomistic and Discrete Dislocation plasticity Plastic deformation and fracture of ductile materials involves important physical phenomena at multiple length scales. Some phenomena (e.g., dislocation nucleation, mobility, crossslip, crack formation, and growth) are intrinsically atomistic. Atomistic studies can address these unit processes involving a few defects but are usually unable to address largerscale deformation except with supercomputers. Plastic deformation at the micron scale, due to motion and interaction of many dislocations, can be described by representing the dislocation core structure, which is the main idea used behind discrete dislocation dynamic modeling of material. A true coupled multiscale 51 model should be capable of dealing with the necessary atomistic degrees of freedom, including all longrange interactions, within one overall computational scheme. Key and distinguishable features of such a model are that, dislocations exist in the continuum region and that they are mechanically coupled to one another and to the atomistic region. Also such dislocations can be passed between the atomistic and the continuum regions so that the plastic deformation is not confined to the atomistic region. Passing of dislocations requires two main developments: (1) detection of dislocations near the atomic/continuum interface that are selected for passing and (2) a procedure for moving such dislocations across the interface. Shilkrot et al. developed such a method, which they call it as: the coupled atomistic and discrete dislocation (CADD) method. In order to detect the dislocations and pass them to the continuum region, they introduce a “detection band” of triangular elements between atoms, with each element sitting on three different slip planes, as shown in Fig. 12. 52 Figure 12. (a) CADD atomistic/continuum interface (c) One detection element, indicating the 3 slip planes on which it lies (Curtin and Miller 2003) If the dislocation passes along one of these planes, it generates a Lagrangian finite strain in the element, element, which can be used to identify the Burgers vector and slip plane of the nucleated defect. For a given crystal structure and orientation, all the possible slip systems and associated strains are known. The detection algorithm monitors the strains in the detection band elements and 53 compares them to the known dislocation slip strains after each energy minimization step. When the strain in an element corresponds to a particular dislocation, the dislocation core is assumed to reside at the centroid of that element. Use of lagrangian finite strain tensor allows for large deformation and rotations. Because multiple dislocations can pass through one element, it is necessary to keep track of all previous slip activity and consider only the displacements due to new defects. After this recognition step the detection establishes the location of the core, the slip plane, and the Burgers vector of the continuum entity. Then the core is artificially shifted along its slip plane from its location in the detection band to a location across the interface in the continuum region by adding the displacement fields associated with a dislocation dipole. These displacements cancel the original core in the atomistic region and add a new core in the continuum region. Once this core is in the continuum region, it is added to the array of discrete dislocation with care being taken to define the branch cut in the continuum displacement field such that it matches the slip plane along which the dislocation moves. CADD is well suited to problems including nanoindentation, fracture, grain boundary sliding, slip transmission across grain boundaries and void nucleation and growth in crystals. So far, this approach has only been implemented in 2D systems, and it has been shown to agree well with the 2D atomistic calculations for Al. Figure 13(a) shows the sample geometry, the rigid indenter and the fully atomistic region as a small box located near the indenter corner, where dislocation nucleation will occur. 54 Figure 13. (a) Nanoindentation simulation using CADD (b) Predicted forcedisplacement curve (c) Positions of nucleated dislocations (Curtin and Miller 2003) The CADD problem contains only about 3000 atoms/nodes (2800 atoms and 200 FE nodes) whereas the full atomistic problem contains about 100 000 atoms. Under indentation, dislocations are nucleated, move away from the indenter into 55 the continuum, and exert a back stress on the atoms near the indenter corner. The predicted forcedisplacement behavior and the positions of nucleated dislocations are shown in figures 13(b) and (c), and compare favorably to the atomistic results and are also presented. The differences in dislocation position are largely associated with the very shallow energy minima for the dislocations in the continuum region, which leads to relatively large differences in position caused by small errors in the predicted stress fields. This approach was also applied to study growth of an atomic scale void in FCC Ni, due to plastic deformation. Recently there have been significant modifications in the nature of atomistic/continuum coupling at the interface (Shilkrot et al. 2004) due to original method showing spurious effects. Compared to CADD the QC method requires full atomic resolution wherever the CADD discrete dislocations are present. 2.5 Conclusion Currently there appears to be a great need to develop general mathematical and computational methods for a truly seamless multiscale approach for computer simulations of nano and Microsystems. In these early stages of development, the computational techniques are developed within a specified range of space and time scales. It is understood that the transition between one spacetime range to another is carried out by a process of handshaking, that is the information gained from a lower scale is summarized into a finite set of parameters, and passed on to the higher scale. This is 56 acceptable as long as the parameters represent a rigorous reduction of the enormous DOF of a lower length scale into a few generalized DOF. However, the theoretical foundations and computational implementation of a more rigorous process remains unresolved. The smooth transitions between various length scales can be ensured is there could be generalizations of the concept of DOF from those associated with spacetime (e.g. geometry), to those representative of statistical configurations (e.g. conductivity, mobility, etc.). 57 CHAPTER 3 PROBLEM STATEMENT A new computational method, namely, the material point method (MPM) was developed by Sulsky and others (1995,1996) at the University of New Mexico with support from the Sandia National Laboratories. Compared to FEM, MPM has the following advantages: (1) MPM is able to handle large deformation in a more natural manner so that mesh lockup present in FEM is avoided; (2) MPM can easily couple with molecular dynamics (MD) simulations because of the use of material points (similar to atoms used in MD) instead of arbitrary sized elements in FEM; (3) Parallel computation is more straightforward in MPM because of the use of a grid structure that is consistent with parallel computing grids; and (4) Use of the background grid in MPM enables structured adaptive refinement for local interested region. As a result, MPM is a logical choice for coupling MD to continuum for multiscale simulations. However, the conventional MPM discretizes the material using a regular grid with uniform square/rectangular (2D) or cubic/brick (3D) cells. Such a mesh is inefficient in dealing with stress concentration issues and 58 cannot be refined to scale down the size to different levels (Tan and Nairn 2002). Consequently, a hierarchical mesh refinement technique in MPM simulations needs to be developed. This study concentrates on the development of multiscale modeling technique. A seamless MD/MPM coupling approach will be developed in which MD modeling will be used in the local area of interest where the variation is at an atomistic scale and MPM simulation in the rest of the material. In this investigation, multiscale modeling and simulations of material properties will be addressed considering tensile testing. This application is chosen because of the extensive use of this technique for determining a wide range of material properties, including modulus, yield strength, ultimate strength, and strain at fracture. MPM/MD simulations of uniaxial tension will be carried on a silicon specimen at different strain rates to illustrate this methodology. A hierarchical mesh refinement technique will be introduced in MPM to scale down from the continuum level into the atomistic level near the MD/MPM transition region where the MPM points overlap the MD atoms. A seamless coupling is necessary for complete compatibility of both strain and stress for MD and MPM regions. Further, elastoplastic constitutive material model will be implemented in MPM for continuum scale simulation of plasticity material behavior. The framework has been established by Sulsky et al. (1996). Elastic plastic MPM at continuum scale coupled with MD at atomistic scale will aid in simulation of nanocrystalline material behavior. 59 CHAPTER 4 COMPUTATIONAL METHODOLOGY 4.1. Introduction. In coupled simulations of material behavior at different length scales, molecular dynamics (MD) is the preferred technique at atomistic scale while material point method (MPM) was chosen for simulation at continuum scale. A seamless coupling between the two techniques was recognized because of the natural similarity between the two methods. For a given body under analysis, molecular dynamics was used at the localized regions of interest, while the rest of the workmaterial is modeled using continuum mechanics approach. The two length scales are coupled by introducing a handshake region in between them, termed as ‘transition region’. A methodology to scale down the locally interested continuum region, to handshake with the atoms in the transition region, was achieved using hierarchical mesh refinement technique. The theoretical background along with computational methodologies adopted for coupled simulations are described in the sections that follow. 60 4.2 Molecular Dynamics (MD) 4.2.1 Theoretical Background In molecular dynamics (MD), the Newtonian second law for the time evolution of a system is expressed by (refer Komanduri et al. 2001,2003) 2 2 mid ri / dt fi , (i = 1,…, N) (4.1) where N is the total number of atoms, mi is the mass of the ith atom with position vector ri, i i 1 2 N f V(r ,r ,...,r ) is the force acting on the ith atom due to interaction with other atoms in the system, and 1 2 N V(r ,r ,...,r ) is the interaction potential. The initial positions and velocities of the atoms, together with the interaction potential, define the whole set of thermodynamic, elastic and mechanical properties of the material. Set of 3N secondorder differential equations [Eq. (4.1)] is often solved by recasting it as a set of 6N firstorder Hamiltonian equations of motion, thus i i i 1 2 N dp / dt f V(r ,r ,...,r ) , i i i dr / dt p /m . (4.2) where pi is momenta. Under the given initial positions and momenta of the system, integration of Eq. (4.2) yields the total trajectory of the system. With the knowledge of the trajectories of all the atoms, one can calculate spatial and temporal distributions of energy, temperature, and pressure, as well as monitor the structural and phase changes in the system. 61 The atomistic potential for the interaction between SiSi atoms herein is a general form of Tersofftype threebody potential (Tersoff 1988,1989). The total energy as a function of atomic coordinates is given by i j ij W 2 1 E , (4.3) where i and j label the atoms of the system and ij W is the bond energy of all the atomic bonds in the control volume represented by W f (r )[f (r ) b f (r )] ij c ij R ij ij A ij , (4.4) where ij r is the interatomic distance. The R f and A f functions are the repulsive and attractive forces between atoms i and j, respectively. These functions are expressed by f (r ) A exp( r ) R ij ij ij ij , (4.5) and f (r ) B exp( r ) A ij ij ij ij . (4.6) The f (r ) c ij function acts as a cutoff function that smoothly attenuates the interaction as ij r approaches the cutoff radius ij S . ij ij ij ij ij ij ij ij ij ij ij c ij 0 for r S ) for R r S S R r R cos( 2 1 2 1 1 for r R f (r ) (4.7) The ij b function is given by i i 2ni 1 n ij n ij ij i b (1 ) , (4.8) where 62 k i, j ij c ik ijk f (r )g( ) , (4.9) with 2 i ijk 2 i 2 i 2i 2i ijk d (h cos ) c d c g( ) 1 , (4.10) where ijk is the bond angle between bonds ij and ik. The parameters used in the Tersoff’s potentials for silicon, namely, A, B, R, S, , , , n, c, d and h are given in table 1 (Tersoff 1988,1989). Table 1. Material parameters in Tersoff’s potential for silicon (Tersoff 1988,1989) Parameter A (eV) B (eV) (Å1) (Å1) n c d h R (Å) S (Å) Silicon 1.83 103 4.71 102 2.48 1.73 1.1 105 0.79 1.0 105 16.22 0.6 2.7 3.0 4.2.2. Strain and Stress In MD simulations, stress and strain relation is usually determined from tensile models (Frankland et al. 2003, Chang and Fang 2003). In each increment of applied deformation, a nominal strain can be determined by the nominal length change normalized by initial length for the entire MD model in the loading direction. In thermodynamics, the stress in a solid is defined as the partial derivative of the internal energy with respect to the strain per unit volume. In continuum mechanics, the components of the stress tensor for a hyperelastic material are given by 63 ij S ij E V 1 , (4.11) where V is the volume of the solid, E is the total internal energy, ij is the components of the strain tensor, and S is entropy. In special cases, if the internal energy is the strain energy, the Hooke’s law is the result of Eq. (4.11). In molecular dynamics, the total energy is the summation of the energy of individual atoms, E , and is expressed as 1 2 E T U M ( ) ( ) 2 v r , (4.12) where T is the kinetic energy,U the potential energy, M the mass, v the magnitude of its velocity, and (r) the potential energy at the atom location r. Thus the stress components for a given atom are written as (M v v F r ) V 1 ij i j i j , (4.13) where V is the atomic volume of atom , i v and j v are components of velocity. i F are the components of force between atoms and from the derivative of the potential energy, and j r are the components of the distance between atoms and . The nominal stress is computed as the average of the atomic stresses for the volume of the model. Thus, the nominal stress components of each model are expressed by (M v v F r ) V 1 ij i j i j , (4.14) where V is the volume of the MD model and equals to V . 64 4.3. The Material Point Method (MPM) 4.3.1 Theoretical Background In the MPM (Sulsky et al. 1995, 1996) a material continuum is discretized into a finite collection of N material points, as shown in Fig. 14. Each material point is assigned a mass (mp, p=1, …, N) consistent with the material density and volume of the point, and all other variables, such as position, acceleration, velocity, strain, and stress. The motion equations are solved at the points for each time step of the analysis using a background computational grid. Information is transferred from the material points to grid nodes by using standard finite element shape functions, n i i i 1 ( ) N ( ) ! x ! x , (4.15) where n is the number of nodes in the grid and subscript i refers to the nodal values of !(x) containing x . Figure 14. MPM grid cells with material points. 65 The Newtonian equations of motion for material points are given by "d2u / dt2 " s "fb , (4.16) where the mass density N p p p 1 m ( ) " # x x , p m is material point mass, p x is material point location, #(xxp) is unit step function, u is the displacement, the acceleration a d2u / dt2 , b f is the specific body force, and s is a specific stress tensor, which is the stress tensor divided by the mass density. The material constitutive relation can be represented as s C(s) and 1 T ( ) ( ) 2 $& u u %' , (4.17) where C(s) is the fourthorder specific elasticity tensor, is the strain tensor. This paper will consider only linear elastic materials at the continuum level. The specific stress components are defined as s / " ij ij (4.18) Thus, Eq. (4.16) can be rewritten as i i s ij, j b a (4.19) Constitutive equations for a linear elastic material are given by k,l s ijkl s ij E v , (4.20a) with E E / " ijkl s ijkl (4.20b) where i v are the velocity components and ijkl E are components of the elasticity tensor. 66 Based on Eqs. (4.16)(4.20), the governing equations in weak form are ( ) () () () " ) " ) " ) T wT dS wb d wa d w d , j s ij i i i (4.21a) and ( E v ) wd 0 k,l s ijkl s ij ( " ) ) (4.21b) where w is the test function which is an arbitrary spatial function. Boundary conditions are u g i on u ) (on displacement boundary) and j s ij i T " n on T ) (on traction boundary), where j n are components of the unit outward normal vector at the boundary surface and u T ) ) ) . It is assumed that w 0 on u ) . A material continuum is divided into a finite collection of discrete infinitesimal regions p ) ( p p 1,...,N ) called material points. Each material point is assigned a mass p m in p ) , where () " ) p m (x)d p and p ) ) . Mass density can then be approximated as a sum of material point masses using a dirac function " # Np p 1 p p (x) m (x x ) (4.22) Applying Eq. (4.22) in Eq. (4.21a, b) converts integrals into summations of quantities evaluated at material points ( ) T p p Np p 1 (p) , j ) p ( , sp ij N p 1 (p) i (p) p N p 1 ) p (i (p) i p wT dS m w b m w a m w (4.23a) and 67 (p) k,l s ijkl ) p ( , sij E v . (4.23b) All variables represented by !(x) , such as the coordinates, displacements, velocities, and accelerations need to be transferred between grid nodes and material points using the shape functions N(x), ! ! N n 1 (x) (n)N(n) (x) . (4.24) where N is the number of nodes in the background grid and !(n) refers to the nth nodal values of !(x) . The numerical solution will be obtained at discrete set of times, t k , where k 1,...,K. Applying the interpolation scheme in Eq. (4.24) to the weak form of momentum equations in Eq. (4.23a) gives * * Np p 1 (p)(n) , j ) p ( k , sp ij n n' 1 N n 1 ) n ( k ) ' n ( ki k(nn') N n 1 wk(n) m a w m N * * N n 1 N n 1 N p 1 k(p) (p)(n) p i k(n) k(n) i k(n) p N b m w Tˆ w (4.25) Applying the interpolation scheme for acceleration into the weak form of Newton’s equation gives n int ext ij j i i j 1 m a f f , (4.26) where the mass matrix is given by N ij p i p j p p 1 m m N ( )N ( ) x x , (4.27) the internal and the external forces on the node i are 68 N int s i p p i p p 1 f m (x ) N (x ) and N ext i i p b p i p p 1 f ˆ m f (x )N (x ) + , (4.28) where i ˆ+ is the surface traction associated with grid point i. The grid point accelerations obtained from Eq. (4.26) are then used to update the position, velocity, stress, strain, and temperature of the material points. 4.3.2. Hierarchical Mesh Refinement (a) (b) Figure 15. Hierarchical mesh refinement in 2D MPM. The conventional MPM method uses a regular grid cell mesh with equally spaced nodes. Adaptive mesh refinement in areas with stress concentration is not possible using such a mesh. To accommodate multiple level meshes, a hierarchical mesh refinement technique is presented herein. For twodimensional 1 2 4 3 5 ,  1 2 4 3 ,  1 2 3 4 5 6 7 8 9 10 Level I Level II Level III 69 situation, the hierarchical mesh refinement with three levels using fournode regular cells (1, 4, 710) and fivenode transition cells (2, 3, 5, 6) are shown in Fig. 15 (a). For the regular cell with four nodes in Fig. 15 (b), the shape functions used in Eq. (4.15) are expressed as (Zienkiewicz et al. 1989,1991), (1 )(1 ) 4 1 N1 ,  , (1 )(1 ) 4 1 N2 ,  , (1 )(1 ) 4 1 N3 ,  , (1 )(1 ) 4 1 N4 ,  . (4.21) where (,, ) are natural coordinates, as shown in Fig. 15 (b). Shape functions for transition cells with five nodes are given by 1 5 N 2 1 (1 )(1 ) 4 1 N ,  , 2 5 N 2 1 (1 )(1 ) 4 1 N ,  , (1 )(1 ) 4 1 N3 ,  , (1 )(1 ) 4 1 N4 ,  , (1 )(1 ) 2 1 N 2 5 ,  . (4.22) Similarly, for threedimensional situation, a hierarchical refinement mesh with two levels is shown in Fig. 16 (a) using regular cells with eight nodes and transition cells with nine (cells 2, 4 and 6) and thirteen nodes (cells 1, 3 and 5). The shape functions for regular cells with eight nodes in Fig. 16 (b) are given by (Zienkiewicz et al. 1989, 1991) (1 )(1 )(1 ) 8 1 N i 0 0 0 ,  , 0 i , ,, , 0 i   , 0 i (i = 1,…,8), (4.23) i( , , ) i i i ,  : 1 (1, 1, 1); 2 (1, 1, 1); 3 (1, 1, 1); 4(1, 1, 1); 5(1, 1, 1); 6(1, 1, 1); 7(1, 1, 1); 8(1, 1, 1), and (,, , ) are natural coordinates as shown in Fig. 16. The shape functions for transition cells from nine to thirteen nodes in Fig. 16 (b) are expressed by 70 (1 )(1 )(1 ) 8 1 Nˆ i i i i ,,  (i = 1,…,8), 1 1 Nˆ N , 2 2 10 11 13 N 4 1 N 2 1 N 2 1 Nˆ N , 3 3 11 12 13 N 4 1 N 2 1 N 2 1 Nˆ N , 4 4 Nˆ N , 5 5 Nˆ N , 6 6 9 10 13 N 4 1 N 2 1 N 2 1 Nˆ N , 7 7 9 12 13 N 4 1 N 2 1 N 2 1 Nˆ N , 8 8 Nˆ N , (1 , )(1 )(1 ). / 4 1 N 2 9 , (1 ,)(1 )(1 ). ,/ 4 1 N 2 10 , (1 , )(1 )(1 ). / 4 1 N 2 11 , (1 ,)(1 )(1 ).,/ 4 1 N 2 12 , (1 , )(1 ).1 / 2 1 N 2 2 13 . (4.24) In 3D MPM computational formulation using the hierarchical mesh with multiple levels, Eq. (4.23) is used for regular cells at each individual level while Eq. (4.24) is employed for the transition cells between neighboring two levels. 71 (a) Regular cell Transition cell (b) Figure 16. 3D MPM regular and transition cells. 1 2 3 5 6 7 11 9 10 12 13  , 4 8 1 2 3 5 6 7  , 4 8 Level I I 6 5 4 Level I 1 2 3 72 4.4. Seamless Coupling between MD/MPM Simulation of dynamic problems using both MD and MPM requires integration of the equations of motion [Eqs. (4.1) and (4.16)]. The format of these two equations is similar, leading to a natural coupling of the two types of simulations. There are essentially three regions in coupled MD/MPM simulations. They are MD, MPM, and their transition region. The MD region is modeled using the approach described in Section 4.2. In MD region, atoms are initially placed at diamond lattice points. The MPM region is modeled using the method described in Section 4.3. In the MPM region, eight material points are equally positioned in each cubic cell. The coupling of the two descriptions of the media is brought about using a transition region where both MPM points and MD atoms are initially overlapped and coincide with the atomic lattice. The width of the transition region is equal to the cutoff distance of the interaction potential used in the MD region. This provides a complete set of neighbors within the interaction range for all atoms in the MD region. Atoms/points that belong to the transition region not only interact via the interaction potential with the MD region but also are part of points in the MPM region. The positions and velocities of these atoms/points in the transition region (both on the MD side and on the MPM side) must be consistent with each other. The atom velocities in the transition region due to the interaction with the MD region furnish the loading boundary condition on overlapped points for the MPM simulation. To avoid a density mismatch at the boundary, the mass at each point for MPM in the transition region is set equal to the mass of each MD atom. 73 Using the above method to deal with the transition region, the displacement as well as strain fields are compatible, while stress fields possibly are not compatible. To ensure complete compatibility of both strain and stress for seamless coupling between MD and MPM, we use the nonlocal elasticity theory (Abraham et al. 1998) in this study. The elastic energy is expanded in terms of the strain in the Taylor series form as [ E /( )] ... 6 1 [ E/( )] 2 1 E( ) E(0) ( E / ) ij kl mn 0 ij kl mn 3 ij kl 0 ij kl 2 ij 0 ij (4.25) Table 2. Material properties of single crystal silicon (Runyan 1999) Material properties Lattice constant (Å) Density (g cm3) Young’s modulus (GPa) C11 (GPa) C12 (GPa) C44 (GPa) Poisson’s ratio Bulk modulus (GPa) Fracture stress (MPa) Silicon 5.43 2.33 130 165.78 63.94 79.57 0.28 100 62 Therefore, the material elastic constants are required to satisfy the following conditions: ij ij 0 c ( E / ) , (firstorder elastic constants) ij kl 0 2 ijkl c [ E /( )] , (secondorder elastic constants) ij kl mn 0 3 ijklmn c [ E /( )] . (thirdorder elastic constants) 74 Using the natural state as the reference configuration, the firstorder elastic constants are zero and only the second and higherorder elastic constants need to be assigned to the continuum. As linear elasticity theory is used in this investigation, the secondorder elastic constants of the lattice and continuum regions require matching. Material properties of single crystal silicon are listed in table 2 (Runyan 1999). 75 CHAPTER 5 NUMERICAL IMPLEMENTATION OF SEAMLESS COUPLING BETWEEN MPM AND MD In this work, a tension model is used to conduct coupled MD/MPM simulations. Using the seamless coupling algorithms described in Chapter 4, a coupled MD/MPM code has been developed for multiscale simulations. The computational scheme or the flow chart used for the coupled MD and MPM simulation of silicon under tension is shown in Fig. 17. To generate the tension model, initial position data need to be assigned first to all atoms and points in both MD and MPM regions, which include atoms/points in the MD/MPM transition regions. Loading and boundary conditions are also prescribed for the entire system. Based on the given temperature, MD computation is first initiated so that the initial velocities for all atoms can be determined. The initial velocities will provide boundary conditions for the MPM region through the transition region. Then, MPM computation can be started and MPM results in the transition region provide boundary conditions for the MD region. Thus, MD computation can be iterated and MD results in the transition region provide new boundary conditions in MPM computation for the next time step. This iteration process continues until the completion of the simulation. For the purpose of illustrating the method, a 76 schematic of simple tension model, which consists of one MD region, two transition regions, and two MPM regions with three hierarchical mesh refinement levels, is shown in Fig. 18. Figure 17. Coupled MD and MPM computational flowchart. ti+1 = ti + 0t ( i = 0, 1, 2,….., n1) if ti+1 = tn Input data 1. Initial position r0 for atoms in MD region, atoms/points in transition regions, and points in MPM regions 2. Loading and boundary condition for entire system v 3. Temperature T for MD. Initial velocities vi0 for atoms in MD region and transition regions t = ti ( i = 0, 1, 2,….., n), 0t MPM computation ri+1, vi+1 for points/atoms in transition regions ri+1, vi+1 for points in MPM regions MD computation ri+1, vi+1 for atoms/points in transition regions ri+1, vi+1 for atoms in MD region END 77 Figure 18. Schematic of tensile model for coupled MD/MPM simulations. The next item of interest is coupling details in the transition region. For the MPM cells in the MD/MPM transition regions (shown in Fig. 18), MPM points are arranged following MD atoms at the atomic lattice sites, which have the interaction with the MD region. In other words, MD atoms and MPM points in the transition region are in onetoone correspondence and are initially overlapped at the same locations, as shown in Fig. 19 (a). It may be noted that MD time step increment is normally much smaller than that in MPM, so that MD and MPM computations have to be matched well in coupling simulations. We use adaptive time step increments in the transition region wherein the MD time step increment and MPM time step increment are the same. For continuum region, the time step increment gradually increases for the area farther away from the atomic region. In this simple model, we chose the same time step increment for both MD and MPM computations. For time integration, the time step increment is 0t = 1.018 fs MD/MPM transition region Atoms Points MPM MD region MPM 78 (1015 s). Details of the steps involved in coupled MD/MPM simulations are given below: (a) (b) (c) Figure 19 (a). MD atoms and MPM points initially overlapped in the transition region. (b). MPM points moving under velocity. (c). MD atoms moving same positions as MPM points. Black: MD atoms Grey: MPM points Dash line: MPM grid Step 1: Initial positions 0 r for all atoms in the MD region and all points in the MPM regions, which include atoms/points in the MD/MPM transition regions, are given according to the dimensions of the tension model and the crystal structure of the material. At the absolute temperature T, MD computation is initiated and initial velocities i0 v for all MD atoms are calculated. The initial velocities at atoms are applied to the overlapped MPM points in the transition regions as a boundary condition for MPM regions, as shown in Fig. 19 (b). Step 2: MPM computation is conducted at the time step t0 and new positions 1 r and velocities 1 v are obtained for all MPM points. The new positions r0 vi0 r1 r1 r0 v1’ 79 at points are applied to the overlapped MD atoms in the transition regions as a boundary condition for MD simulation as shown in Fig. 19 (c). Step 3: MD computation is performed at the time step 0 t and the updated positions r1 ' and velocities 1 v ' are obtained for all atoms. The updated velocities at atoms are applied to the overlapped MPM points in the transition regions as a new boundary condition for MPM simulations at the next time step 1 t . Step 4: Repeat the same algorithm to steps 2 and 3 in turns until the entire computation for the tension model is implemented by means of coupled MD and MPM approach. During coupling of MD with MPM simulations, MPM always provides the updated positions for MD atoms in the transition regions, while MD in turn provides the updated velocities for MPM points in the transition regions. It should be noted that following this algorithm from step 1 to step 4, both atoms for MD and points for MPM in the transition region naturally maintain displacement compatibility so that strains are compatible as well. However, the stress compatibility cannot be ensured automatically. Additionally, mass mismatch exists for both MD atoms and MPM points in the transition region, since mass determined by the mass density and the occupied volume is assigned to each point in the MPM region. To ensure complete compatibility of both strain and stress for seamless coupling between MD and MPM, nonlocal elasticity theory is 80 used in this study. Based on this theory, for homogenous, linear elastic materials, only the secondorder elastic constants require matching between atomistic and continuum regions. In this investigation, we take elastic constants for continuum region from the lattice theory and use the same mass for both MD atoms and MPM points in the transition regions. 81 CHAPTER 6 MULTISCALE SIMULATIONS OF TENSILE TESTING 6.1. Introduction In this investigation, multiscale modeling and simulations of material properties are addressed considering tensile testing. This application is chosen because of the extensive use of this technique for determining a wide range of material properties, including modulus, yield strength, ultimate strength, and strain at fracture. In this work, a tension model is used to conduct coupled MD/MPM simulations. Using the seamless coupling approach, multiscale simulations of tensile testing are carried on single crystal silicon workmaterial. Results are presented in terms of stress  strain relationships at several strain rates, as well as the rate dependence of uniaxial material properties. The simulation details and the results obtained are given in the sections that follow. 82 6.2. Coupled MD/MPM tensile testing model. Fig. 20 is 3D animation of coupled MD/MPM model at the initial stage of tensile testing showing the MD region, MD/MPM transition regions, and MPM regions. There are three levels in the MPM region. The dimensions of this model are 32.6 Å 32.6 Å 222.8 Å. There are 2,901 atoms in the MD region and 9,280 material points in the MPM regions, which include atoms/points in the MD/MPM transition regions. In the transition regions, MD atoms and MPM points are initially overlapped in onetoone correspondence and coincide with the atomic lattice. In the MD region, the crystal is set up with a cube orientation [001] and the tensile loading is applied along the [001] direction under controlled velocity. The initial temperature is 293 K for the MD simulations. The elastic constants, the Young’s modulus, the Poisson’s ratio, and other relevant parameters, given in Table 2, are used in the MPM regions. The coupling simulations on the tension model are conducted for the silicon workmaterial at different strain rates, namely, 44.88 ns1, 53.86 ns1, 62.84 ns1, 67.33 ns1, 80.79 ns1, and 89.76 ns1. It may be noted that very high strain rates are used of necessity to minimize the computational time for coupling simulations. For the coupling simulation at a strain rate of 44.88 ns1, the total simulation time was 4 ps (1012 s) and the time step increment used was 1.018 fs. It took ~ 3 hrs of computational time using a PC with a 2.8 GHz processor. 83 Figure 20. Animation of coupled MD/MPM model showing initial stage of tensile testing. MD Atoms MPM Points MPM Points MD region Level I MD/MPM transition region Level II Level III MPM region Level III Level II Level I MPM region 84 6.3. Deformation and failure mechanism under tensile loading Fig. 21 shows snapshots of coupling simulation for the entire model at various stages of tensile testing at a strain rate of 44.88 ns1 and Fig. 22 shows closeup views of Fig. 21. The animations of coupled MD/MPM simulation show that the tensile specimen experiences with increasing elongation in elasticity, generation of dislocations and plasticity by slip, void formation and propagation, formation of amorphous structure, necking, and final fracture in the MD region while only elastic deformation is present in the MPM regions at either end. In the initial stage, there is uniform stretching with elastic deformation and no dislocations occurring in the crystal lattice, as shown in Fig. 23 (b). The limit of elasticity is reached when the strain reaches approximately 0.1, as shown in Fig. 27. With further increase in deformation, the dislocations are seen to be generated locally with their densities gradually increasing in the plastic deformation region. Slip takes place along [111] at about 45o to the loading direction as shown in Fig. 23 (c). Plasticity by slip lasts until the formation of voids and semiamorphous structure, as shown in Fig. 23 (d). Ultimate tensile stress is reached at a strain of ~ 0.3, as shown in Fig. 27. Further increase in elongation results in local interactions or rearrangements of silicon atoms, resulting in somewhat amorphous structure [see Figs. 23 (e) and (f)]. The deformation is then concentrated in a narrow region resulting in necking, as shown in Fig. 21 (d). Finally, separation of the tensile specimen occurs as shown in Fig. 21 (e). 85 (a) (b) (c) (d) (e) Figure 21. Snapshots of MD/MPM coupling simulation at strain rate 44.88 ns1. 86 ( a ) ( b ) ( c ) ( d ) ( e ) Figure 22. Snapshots of MD/MPM coupling simulation at strain rate 44.88 ns1 (magnified view). 87 Figs. 22 (a) to (e) are 3D animations of MD simulation at higher magnification (of Fig. 21) showing different stages of tensile testing at a strain rate of 44.88 ns1. Slip taking place at 45o to the loading direction can be clearly seen [Figs. 22 (b) to (c)]. The formation of amorphous structure in the central region and void formation can also be clearly seen in Fig. 22 (d). Fig. 22 (e) shows concentration of amorphous structure in the central region followed by necking and separation. This corresponds to the steepest slope in Fig. 27 leading to final rupture of the sample. Figs. 23 to 27 show animations in the midsection at different stages of tensile testing at different strain rates from 44.88 ns1 to 89.76 ns1. Initially the deformation is predominantly elastic. This corresponds to steep elastic slope in Fig. 27. This is followed by slip on the preferred [111] direction. Further increase in strain leads to the formation of inhomogeneous structure leading to amorphization in the entire necked region [see atoms marked by circles in Fig. 23 (d)]. A similar amorphization can be seen in Figs. 24 (d), 25 (d) and 26 (d). Ultimately, separation of the sample takes place with drastic reduction in stress. With increase in strain rate, the crystallinity of the silicon is maintained up to larger strains. Also, with increase in strain rate, the strain at fracture increases. It may be noted that necking that occurs before final failure for silicon is significantly less as compared to ductile materials, such as copper, aluminum (Komanduri et al. 2001). The catastrophic failure leads to a sudden drop in the stressstrain curves as shown in Fig. 27. 88 (a) (b) (c) (d) (e) (f) Figure 23. Snapshots of MD/MPM coupling simulation at strain rate 44.88 ns1. 89 (a) (b) (c) (d) (e) (f) Figure 24. Snapshots of MD/MPM coupling simulation at strain rate 67.33 ns1. 90 (a) (b) (c) (d) (e) (f) Figure 25. Snapshots of MD/MPM coupling simulation at strain rate 80.79 ns1. 91 (a) (b) (c) (d) (e) (f) Figure 26. Snapshots of MD/MPM coupling simulation at strain rate 89.76 ns1. 92 6.4 Effects of strain rate on material properties The effects of strain rate on material properties have also been investigated for the coupled MD/MPM simulations. Under different strain rates, slip essentially takes place at about 450 with respect to the loading direction followed by the formation of amorphous structure and generation of voids in the necked region and their coalescence, leading to final failure. The engineering stressstrain curves at different strain rates are shown in Fig. 27. They show that at different strain rates material behaves differently at different stages including elastic and plastic deformation, and final rupture. Thus, the strain rate affects material properties significantly. However, before the strain reaches ~ 0.1, there is no dislocation generation and the initial part of engineering stressstrain curve is essentially elastic at different strain rates in the range from 44.88 ns1 to 89.76 ns1. From the slope of the stressstrain curves, the Young’s modulus is determined as ~ 121 GPa, which is very close to that obtained from the lattice theory (Runyan, 1999). As expected, the elastic modulus is nearly constant and is not affected by the strain rate. The curve in Fig. 28 (a) shows the effect of the strain rate on the tensile strength for silicon. The tensile strength increases from ~ 15 to ~ 18 GPa with increase in the strain rate from 44.88 ns1 to 89.76 ns1. Fig. 28 (b) shows the variation of the ultimate strain with strain rate. It shows that the ultimate strain increases rapidly up to ~ 0.54 as the strain rate increases from 44.88 ns1 to 67.33 ns1, and then increases gradually up to ~ 0.57 as the strain rate increases from 67.33 ns1 to 89.76 ns1. The ultimate strain vs. the strain rate curve 93 (GPa) 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0 5 10 15 20 53.86 ns 67.33 ns1 1 44.88 ns1 62.84 ns1 80.79 ns1 89.76 ns1 Figure 27. Engineering stressstrain curves from coupling MD/MPM simulations for silicon. changes slope at a strain rate of ~ 67.33 ns1. Observation of the animations indicates that defect initiation and coalescence are somewhat different across the strain rate. At lower strain rates, i.e. below 67.33 ns1, [for example at 44.88 ns1, as shown in Fig. 23(d)], voids are formed in the central region of the tensile specimen. These voids coalesce to larger defects resulting in necking and ultimate failure. While at higher strain rates (67.33 ns1 or higher), [for example at 89.76 ns1, as shown in Figs. 26(d) and (e)], defects are formed near the surface at about ¼ and ¾ of the MD region (see regions marked by circles in Fig. 26). These defects grow in size and cover the entire area, leading to decreased load bearing capacities at these locations, resulting in some bulging, as shown in Figs. 26 (d)  (e). Compared with failure at the lower strain rates 94 40 50 60 70 80 90 12 14 16 18 20 (ns ) 1 u (GPa) (a) Tensile strength vs. strain rate. 40 50 60 70 80 90 0.3 0.4 0.5 0.6 0.7 (ns ) 1 u (b) Ultimate strain vs. strain rate. Figure 28. Effects of strain rate on material properties. 95 [e.g., Fig. 23(d)], there are more defect sites at higher strain rate [Fig. 26(e)]. Once the defects at multiple sites coalesce, necking occurs, leading to a quick rupture of the material. The change in the slope of the ultimate failure strain shown in Fig. 28 (b) is most likely the result of the change in defect formation mechanism from relatively low to high strain rates. However, it is not clear at this stage the underlying reason for a change in the defect formation mechanism across the strain rates. Further investigation is needed to explore this feature. 96 CHAPTER 7 ELASTOPLASTIC ANALYSIS USING MATERIAL POINT METHOD 7.1.1 Metal plasticity models In processing of metallic materials, elasticplastic constitutive material behavior is involved. Numerical analysis for this class of materials is an important aspect. A wide range of processes, such as metal cutting, metal forming involve large deformation and plastic flow of the material. The material point method (MPM) has been recently found useful for such applications (Sulsky et al. 1994). In plasticity, numerous constitutive models exist and the main choices are between rate independent and rate dependent plasticity, between the von Mises and Hill’s yield surface. Further, we have choice within rate independent model as either isotropic hardening or kinematic hardening. Rate independent model is adequate for modeling the response of metals at low strain rates, and at operating temperatures much lower than the melting temperature. For metals, yielding is independent of the hydrostatic stress. Experiments advocate that von Mises yield criterion provides somewhat better fit to experimental data than Tresca yield criterion. Hence von Mises yield surface is used in this formulation to define isotropic yielding. A state of stress on the von 97 Mises yield surface causes plastic deformation and yielding. In case strain hardening occurs, the radius of the yield surface increases. In the case where polycrystalline metals obey isotropic plastic hardening law, the radius of the yield surface changes in all directions. In dynamic problems such as tension and nanoindentation, which involve finite strains, as well as in manufacturing processes, such as metal forming and metal cutting, large plastic strains are present and plastic strain does not continuously reverse direction sharply. In such applications, isotropic hardening theory is usually considered. Associated plastic flow rule is used, which means that the plastic deformation increment has the same direction as the outward normal to the yield surface. It should be noted that, there are also some limitations posed by the model, namely, the model does not predict behaviour under cyclic loading correctly. Since the incremental plasticity model is based on the assumption of small displacements and small strains, it will not predict plastic strains correctly if the principal axes of stress rotate significantly during deformation. The general equations for the plasticity model using von Mises yield surface, isotropic hardening and associative flow rule, are summarized in the following section. Two examples, one, a three  dimensional tensile bar under impact load, and the other Taylor Impact bar problem are used to validate the code. For the later, results are presented to demonstrate axisymmetric formulation of the material point method (Sulsky and Schreyer, 1996). In computational mechanics finite element method (FEM) is being used for a wide range of problems with satisfactory results. Realizing the similarity in 98 basic dynamic equations that are solved by both MPM and FEM and that too using explicit formulation, FEM results are chosen for comparison and validation of MPM elastoplastic constitutive material model. 7.1.2 Summary of rate independent multiaxial plasticity model with von Mises yield surface and isotropic hardening. 1 Using Incremental approach and additive strain decomposition: pi j e ij ij , dt d dt d E dt d plm mn n ijmn ij , where, ij : Total incremental 2nd order strain tensor. t : time. e ij : Elastic part of the ij . pi j : Plastic part of the ij . Eijmn: 4th order elasticity tensor. ij : 2nd order stress tensor. 1 von Mises yield criterion, f( ij ) : 0 p ij SijSij H 2 3 f( ) , ij ij kk ij 3 1 S , where, ij S : Deviotoric part of ij . kk : 1st invariant of ij . 99 p : Equivalent plastic strain. p H : Radius of the yield surface. ij : 2nd order unit tensor. Also, pi j pi j p , where, pi j : Components of plastic strain tensor. 1 Isotropic strain hardening law: Linear hardening: . p / H = 0 + c p . Power law hardening: . p / H = 0 + c . /p m . where, 0 : Virgin yield stress in uniaxial stressstrain behavior. m: material constant. c: Slope of yield stress versus plastic strain curve, positive and negative values of which corresponds to hardening or softening, respectively, of the material. 1 Associative plastic flow rule: mn mn ij ij pi j S S S dt d d df dt d dt d 2 2 , where, dt: time increment. 2 : Measure of the plastic strain path length, and is a nonnegative parameter. 1 Loading / unloading condition in Kuhn  Tucker complimentary form: Find trial stress increment, 100 0 trial = E 0 and initialize 0 0 . If f( ) 0 and 0 d df 0 3' % 4& $ , then instantaneous elastic process, else If f( ) = 0 and 0 5 0 3' % 4& $ , then instantaneous plastic process. The loading / unloading condition above can be solely determined using trial state only, which is ultimately implemented in return mapping algorithm. 7.1.3 Radial return mapping algorithm for threedimensional nonlinear isotropic hardening plasticity. 7.1.3.1 Main algorithm. 1) Compute trial elastic stress and deviatoric stress, assuming all the strain increment is elastic. 6tr. /7 I 3 1 n 1 n 1 n 1 e , . / pn n 1 trial n 1 s 2 e e , where, I : Unit tensor; n 1 e : deviatoric part of total strain. n 1 : strain at current time; trial1 n s : trial deviatoric stress; K: Material bulk modulus; : shear modulus; p n e : plastic part of total strain. Subscripts (n+1) and (n) indicate current and previous time step values respectively. 101 2) Check yield condition. . p / n trial1 n trial n 1 H 3 2 f s IF f trial 0 n 1 then all the stress increment is elastic so update total stress ( ) and exit. . / trial1 n 1 n 1 n Ktr I s , exit. ELSE plas
Click tabs to swap between content that is broken into logical sections.
Rating  
Title  Multiscale Simulation from Atomistic to Continuum  Coupling Molecular Dynamics (MD) with Material Point Method (MPM) 
Date  20041201 
Author  Daphalapurkar, Nitin Pandurang 
Keywords  Multiscale simulations, atomistic/continuum coupling, molecular dynamics, material point method, tensile testing, single crystal silicon 
Department  Mechanical Engineering 
Document Type  
Full Text Type  Open Access 
Abstract  Failure in single crystals and polycrystalline materials usually involve processes such as dislocation, cleavage, macrocrack initiation and growth as well as coalescence until final fracture. Multiscale modeling is necessary to understand the mechanical behavior of materials from atomistic to continuum scales. MPM has been used for continuum simulation. The use of material points at the continuum level provides a natural connection with the atoms in the lattice at the atomistic scale. A hierarchical mesh refinement technique in MPM is presented to scale down the continuum level to the atomistic level, so that material points at the fine level in MPM are allowed to directly couple with the atoms in the MD. A onetoone correspondence of MD atoms and MPM points is used in the transition region, and nonlocal elastic theory is used to assure compatibility between MD and MPM regions, so that seamless coupling between MD and MPM can be accomplished. A single crystal silicon workmaterial under uniaxial tension is used in demonstrating the viability of the technique. A Tersofftype, threebody potential was used in the MD simulations. Further, elastic plastic constitutive material model is integrated with threedimensional MPM to aid simulation of nanocrystalline material behavior at continuum scale. A new multiscale simulation approach is introduced that couples atomistic scale simulations using MD with continuum scale simulations using MPM. The coupled MD/MPM simulations show that the silicon under nanometric tension experiences with increasing elongation in elasticity, dislocation generation and plasticity by slip, void formation and propagation, formation of amorphous structure, necking, and final rupture. Results are presented in terms of stress  strain relationships at several strain rates, as well as the rate dependence of uniaxial material properties. This new multiscale computational method has potential for use in cases where a detailed atomisticlevel analysis is necessary in localized spatially separated regions whereas continuum mechanics is adequate in rest of the material. 
Note  Thesis 
Rights  © Oklahoma Agricultural and Mechanical Board of Regents 
Transcript  MULTISCALE SIMULATION FROM ATOMISTIC TO CONTINUUM  COUPLING MOLECULAR DYNAMICS (MD) WITH MATERIAL POINT METHOD (MPM) By NITIN P. DAPHALAPURKAR Bachelor of Engineering University of Pune Pune, Maharashtra, India 2002 Submitted to the faculty of the Graduate College of the Oklahoma State University in partial fulfillment of the requirements for the Degree of MASTER OF SCIENCE December 2004 ii MULTISCALE SIMULATION FROM ATOMISTIC TO CONTINUUM  COUPLING MOLECULAR DYNAMICS (MD) WITH MATERIAL POINT METHOD (MPM) Thesis Approved: Dr. Hongbing Lu Dr. Ranga Komanduri Dr. Samit Roy Dr. A. Gordon Emslie Advisor Dean of the Graduate College iii ACKNOWLEDGEMENTS First, I would like to express my sincere gratitude and appreciation to Professor, Dr. Hongbing Lu for his guidance, support, and encouragement during the entire course of this study. Working under his advisement was really a valuable experience for me. I would also like to express my appreciation to Dr. S. Roy and Dr. R. Komanduri for their constant guidance and encouragement throughout this study. Thanks are due to Dr. K. A. Salam for reviewing my thesis work. I would like to express my appreciation to Dr. Bo Wang for his invaluable guidance and direct contribution towards the completion of this thesis. Without their great friendship, understanding, and encouragement, it would be hard to believe that I could come this far. My master studies under them have been a wonderful and unforgettable period. The work was supported by a grant from the Air Force Office of Scientific Research (AFOSR) through a DEPSCoR grant (No. F496200310281). The author thanks Dr. Craig S. Hartley, Program Manager for Metallic Materials Program at AFOSR for his interest and support of this work. Thanks are also due to M. Malshe, C. Nagasubramanian, R. Narulkar, V. Karuppiah, D. Joshi, J. Ma, H. Viswanathan as well as to the Molecular Dynamics and FEM research group under Dr. Komanduri, for their support and helpful conversations. iv I wish to extend my gratitude to R. Raghav and Y. Liu for their encouragement and friendship. Special thanks are due to Dr. Gary Young for his support during this graduate study. Finally, I am indebted to my father Prof. Pandurang B. Daphalapurkar, and my mother Rupali P. Daphalapurkar for their firm belief in education and beloved concern all the time. Their encouragement and support all the time have helped me successfully fulfill this goal. v ABSTRACT A new multiscale simulation approach is introduced that couples atomistic scale simulations using molecular dynamics (MD) with continuum scale simulations using the recently developed material point method (MPM). In MPM, material continuum is represented by a finite collection of material points carrying all relevant physical characteristics, such as mass, acceleration, velocity, strain, and stress. The use of material points at the continuum level provides a natural connection with the atoms in the lattice at the atomistic scale. A hierarchical mesh refinement technique in MPM is presented to scale down the continuum level to the atomistic level, so that material points at the fine level in MPM are allowed to directly couple with the atoms in the MD. A onetoone correspondence of MD atoms and MPM points is used in the transition region, and nonlocal elastic theory is used to assure compatibility between MD and MPM regions, so that seamless coupling between MD and MPM can be accomplished. A single crystal silicon workmaterial under uniaxial tension is used in demonstrating the viability of the technique. A Tersofftype, threebody potential was used in the MD simulations. The coupled MD/MPM simulations show that silicon under nanometric tension experiences with increasing elongation in elasticity, dislocation generation and plasticity by slip, void formation and vi propagation, formation of amorphous structure, necking, and final rupture. Results are presented in terms of stress  strain relationships at several strain rates, as well as the rate dependence of uniaxial material properties. Isotropic elastoplastic constitutive model is integrated with the three dimensional formulation of MPM to simulate the material behavior of nanocrystalline materials at continuum scale. Validation of the constitutive model is done using two examples, impact load on a tensile bar and Taylor impact problem. This new multiscale computational method has potential for use in cases where a detailed atomisticlevel analysis is necessary in localized spatially separated regions whereas continuum mechanics is adequate in rest of the material. vii TABLE OF CONTENTS Chapter Page 1. INTRODUCTION………………………………………………………………… 1 1.1. Need for multiscale modeling…………………………………………….. 1 1.2. Aspects of multiscale simulation…………………………………………. 2 1.3. Overview……………………………………………………………………. 5 2. LITERATURE REVIEW…………….…………………………………………… 7 2.1. Introduction………………………………………………………………… 7 2.2. Coupled atomistic and continuum simulation approach………………. 11 2.3. Multiscale modeling of plasticity using coupled mesoplasticity and continuum simulation……………………………………………………... 38 2.4. Coupled Atomistic and Discrete Dislocation plasticity……………… .. 50 2.5. Conclusion……………………..…………………………………………. 55 3. PROBLEM STATEMENT…………….………………………………………… 57 4. COMPUTATIONAL METHODOLOGY………………………………………... 59 4.1. Introduction………………………………………………………………… 59 4.2. Molecular Dynamics (MD)………………………………………………… 60 4.3. The Material Point Method (MPM).……………….…….............………. 64 4.4. Seamless Coupling between MD/MPM.............................……………. 72 5. NUMERICAL IMPLEMENTATION OF SEAMLESS COUPLING BETWEEN MPM AND MD…………………….……………………………………………. 75 viii 6. MULTISCALE SIMULATIONS OF TENSILE TESTING………………….… 81 6.1. Introduction………………………………………………………………… 81 6.2. Coupled MD/MPM tensile testing model………………………………. 82 6.3. Deformation and failure mechanism under tensile loading…………… 84 6.4. Effects of strain rate on material properties……………………………. 92 7. ELASTOPLASTIC ANALYSIS USING MATERIAL POINT METHOD……. 96 7.1. Metal Plasticity models…………………………………………………… 96 7.2. Numerical results………………………………………………………….. 103 8. CONCLUSIONS AND FUTURE WORK ………………………………………116 BIBLIOGRAPHY......……………………………………………………………..118 ix LIST OF TABLES Table Page 1. Material parameters in Tersoff’s potential for silicon (Tersoff 1988,1989)… 62 2. Material properties of single crystal silicon (Runyan 1999)……………… 73 x LIST OF FIGURES Figure Page 1 (a). Schematic of simulations at various levels (Komanduri, 2002).……… 3 (b). Schematic of various length and time scales in tribology (Hutchings, 2002)…………………………………………………………. 3 2. Representative illustration of FE/MD handshake at the interface (Rudd and Broughton, 2000)…………………………………………………. 13 3. The stress waves propagating through the slab using a finely tuned potential energy color scale (Abraham et al. 2000)………………………………………. 14 4. Transition region for QC method (Curtin and Miller 2003)….……………. 17 5. (a) FE mesh used to model dislocationGB interaction using QC method (Shenoy et al. 1999)………………………………………………………………. 20 (b) Snapshots of atomic positions at different stages in the deformation history. (Shenoy et al. 1999)…………………………………………………. 21 6. Quasicontinuum simulation of fracture at atomic scale (Miller et al. 1998) 23 7. Nanoindentation simulation on thin film Al using QC method (Tadmor et al. 1999)……………………………………………….………… 25 8. Schematic diagram of the construction of CWM (a) from the LJ system of 1x1 units of length discretized into 256x256 points (b) from the Potts system of 2x2 units of length discretized into 128x128 points. (Frantziskonis and Deymier, 2000)…………………………………………. 31 xi 9. Connection between continuum region and atomistic region (Tan 2001).. 36 10. Nanoindentation simulation stress contour using multiscale model of plasticity. (Zbib and Rubia 2002)…………………………………………………. 43 11. Results of nanoindentation simulations using multiscale plasticity model. (a) Developed dislocation structure underneath the indenter. (b) Predicted indentation depth showing the development of surface distortions. (Zbib and Rubia 2002)……………………………………………………… 44 12. (a) CADD atomistic/continuum interface…………………………………… 52 (b) One detection element, indicating the 3 slip planes on which it lies (Curtin and Miller 2003)…………………………………………………… 52 13. (a) Nanoindentation simulation using CADD……………………………… 54 (b) Predicted forcedisplacement curve……………………………………. 54 (c) Positions of nucleated dislocations (Curtin and Miller 2002)………… 54 14. MPM grid cells with material points……………………………………….. 64 15. Hierarchical mesh refinement in 2D MPM……………………………….. 68 16. 3D MPM regular and transition cells……………………………………… 71 17. Coupled MD and MPM computational flowchart………………………… 76 18. Schematic of tensile model for coupled MD/MPM simulations………… 77 19 (a) MD atoms and MPM points initially overlapped in the transition region……………………………………………………………………… 78 (b). MPM points moving under velocity……………………………………….. 78 (c). MD atoms moving same positions as MPM points…………………….. 78 xii Figure Page 20. Animation of coupled MD/MPM model showing initial stage of tensile testing…………………………………………………………………………. 83 21. Snapshots of MD/MPM coupling simulation at strain rate 44.88 ns1…. 85 22. Snapshots of MD/MPM coupling simulation at strain rate 44.88 ns1 (magnified view)……………………………………………………………… 86 23. Snapshots of MD/MPM coupling simulation at strain rate 44.88 ns1….. 88 24. Snapshots of MD/MPM coupling simulation at strain rate 67.33 ns1….. 89 25. Snapshots of MD/MPM coupling simulation at strain rate 80.79 ns1….. 90 26. Snapshots of MD/MPM coupling simulation at strain rate 89.76 ns1…. 91 27. Engineering stressstrain curves from MD/MPM coupling simulations for silicon………………………………………………………………………… 93 28. Effects of strain rate on material properties……………………………… 94 (a). Tensile strength vs. strain rate. (b). Ultimate strain vs. strain rate. 29. 3D bar represented by collection of material points …………………….. 103 30. Applied load history…………………………………………………………. 104 31. von Mises vs time plot for elastic material model……………………….. 105 32. von Mises vs time plot for elasticperfectplastic material model with yield stress of 1.5 Pa…………………………………………………………….. 105 33. von Mises vs time plot for elasticplastic material model with power law hardening …………………………………………………………………… 105 34. von Mises stress contour at time = 10.89 s…………………………… 107 35. Stress XX contour at time = 10.89 s…………………………………… 108 xiii Figure Page 36. Displacement along X contour at time = 10.89 s……………………….. 108 37. Equivalent plastic strain contour at time = 10.89 s…………………….. 108 38. Original shape represented by material points at 0 s………………….. 110 39. Deformed shape represented by material points at 24 s……………… 110 40. Deformed shape at 76 s………………………………………………….. 110 (a). MPM points representing the deformed bar. (b). Deformed FEM mesh representing the deformed bar. 41. Equivalent plastic strain contour at 24 s…………………………………..113 (a) Axisymmetric MPM (b) Axisymmetric FEM 42. 3D Equivalent plastic strain distribution contour at material points………114 (a) At time = 24 s (b) At time = 90 s 43. Time history of the total kinetic energy…………………………………….. 115 44. Time history of the free end displacement …………………………… ….. 115 1 CHAPTER 1 INTRODUCTION 1.1. Need for multiscale simulation Material failure at all length scales experiences such processes as elastic deformation, dislocation generation and their propagation, cleavage, crack initiation and growth, and final rupture. Despite significant developments in materials simulation techniques, the goal of reliably predicting the properties of new materials in advance of material characterization and fabrication has not yet been achieved. It is also not possible to predict the properties of the material at nanoscale knowing the properties of the material at macrolevel or vice versa. This situation exists for several reasons that include a lack of reliable potentials to describe the behavior of the material at the atomistic level, computational limitations, inadequate modeling of the process, absence of scaling laws, and difficulties associated with the experimental measurement of properties even at the microscale, let alone at the nanoscale. Multiscale simulation technique can be applied to reliably predict the properties of new materials in advance of material characterization and fabrication. Scaling laws governing the mechanical behavior of materials from atomistic (nano), via mesoplastic (micro), to continuum (macro) are very 2 important to numerous DoD applications, such as the development of a new class of aircraft engine materials, or new steels for naval battle ships, or new tank armor materials for the army, or numerous microelectromechanical components for a myriad of applications, for the following three reasons: (1) information on the mechanical behavior of materials at nanolevel is not presently available as input to nanotechnology for the manufacturing of nanocomponents or microelectromechanical systems (MEMS). For example, nanostructures may possess unique properties in view of their very high surface to volume ratio, or the nanostructures might be relatively free of defects with strength approaching the theoretical values, (2) applications where two length scales of different orders of magnitude are involved. For example, one is atomistic (nano) and the other mesoplastic (micro) as in nanoindentation, or, one is mesoplastic (micro) and the other continuum (macro) as in conventional indentation, and (3) it may be possible to extend the knowledge accumulated over time on material behavior at the macro (or continuum) level to the atomistic (or nano) level, via mesoplastic (micro) level. 1.2. Aspects of multiscale simulation In recent years, combined atomistic and continuum simulations have received increasing attention due to their potential linkage between structureproperty relationships from nano to macrolevels (King el al. 1995, Phillips 1998, Shenoy et al. 1999, Smirnova et al. 1999, Smith et al. 2001, Ogata et al. 2001, and Shilkrot et al. 2002). Various computational methods were developed to 3 Figure 1 (a). Schematic of simulations at various levels (Komanduri, 2002). Figure 1 (b). Schematic of various length and time scales in tribology (Hutchings, 2002). 4 couple spatial and temporal scales. Fig. 1 (a) is a schematic of the various computer simulation methods used in the development of scaling laws as a function of spatial and temporal variables (Komanduri, 2002). At the subatomic level, we use quantum mechanics or ab initio calculations to develop the potentialenergy functions. With increase in length (or time) scale, we model material behavior using molecular dynamics (MD) and Monte Carlo (MC) simulations, then to micro or mesoplastic, and finally to continuum mechanics, such as finite element method (FEM). Fig. 1 (b) is a schematic showing various spatial and temporal scales for a typical tribological phenomenon for a wide range of applications (Hutchings, 2002). It may be noted at the outset that any scaling law developed would involve both temporal and spatial scales. At the atomistic level, MD is the preferred mode while at the continuum level, FEM is the dominant numerical technique used. Consequently, most published results on combined atomistic and continuum simulations have focused on MD/FEM coupling. At Oklahoma State University considerable work is being conducted in multiscale modeling and simulation for materials processing. The main objective of this ongoing research is to develop a fundamental understanding of indentation and tension tests at various scales and develop a computer simulation code that would bridge the gap from atomistic to continuum, via. mesoplastic or microscopic behavior. It is hoped that the outcome of this work would provide the required inputs to Computer Aided Design (CAD) software so that relevant output parameters can be used by the engineer in the design of 5 components at any level of detail, namely, large macro components, or small microcomponents, or even extremely small nanocomponents. It may be noted that the work described, herein, is part of a larger ongoing research to aid materials modeling and simulation. 1.3. Overview In Chapter 2, major multiscale modeling approaches are reviewed. The major focus of the review is on atomistic/continuum coupling. Material point method (MPM) is introduced and its potential to handle a wide range of problems is reviewed and discussed. In Chapter 4, the theoretical background and computational methodologies used in MD and MPM are reviewed. Based on these two simulation techniques, seamless coupling between these two for ability to perform coupled simulations is recognized. Chapter 5 discusses numerical implementation methodology of seamless coupling between MD/MPM. In Chapter 6, simulation approach, material properties, and results of multiscale simulations of tensile testing are presented. In tensile simulations, the effects of strain rate on the nature of deformation and material properties are investigated based on the simulated stress  strain relations for single crystal silicon. 6 Chapter 7 discusses elastoplasticity, and its computational methodology for numerical implementation in MPM. Validation results of the same are presented. Finally in Chapter 8, conclusions are made out of present investigation and offer some suggestions for future work. 7 CHAPTER 2 LITERATURE REVIEW 2.1. Introduction Continuum approach of modeling the behavior of materials, dominated the research method for numerous decades. The analysis using continuum mechanics approach implicitly assumes that, all the atomic scale dynamics and defect evolutions that determine the mechanical properties are averaged over time and large length scales. In fact, this notion stemmed from the desire to reduce the enormous number of degrees of freedom (DOF) in a given system to as minimum as possible. Thus, most of the relevant information in constitutive equations represents some averaged behavior of many, many atoms. However, the processes such as, dynamic fracture instability (Abraham et al. 1994), intersonic crack motion (Abraham and Gao, 2000), dynamic competition of dislocation emission and cleavage (Tan and Yang, 1994, Zhang and Wang, 1995), atomic inertial effect and chaotic atom motion at crack tip (Tan and Yang 1996), dislocation patterns, bifurcation phenomenon, etc., which are seen in material systems at the nano and micro scales, cannot be described by continuum theory alone. In the abovementioned processes, the length and time scales reach those of typical defects: point defects, dislocations, interfaces, and 8 grain boundaries. Hence, in case of all the abovementioned phenomena that represent material catastrophes, the averaging techniques will not yield the correct information. Nonetheless, the dawn of large scale computing is driving the art and science of modeling material phenomena into a tantalizing new direction. Instead of attempting to reduce the complexity of the material system’s behavior by a process of reduction of its DOF, one is trying to represent large DOFs, and solve for them numerically. The result is a real numerical experiment, whose outcome is not known a priori. Whether that makes sense can only be tested by confronting the outcomes of computer simulations with a limited range of experimental observations. In atomistic simulations under the given initial positions and momentum of the system, integration of Newtonian second law equation yields the total trajectory of the system. With the knowledge of the trajectories of all the atoms, one can calculate spatial and temporal distributions of energy, temperature, and pressure, as well as map out all atomistic scale processes. Advances have been made in atomistic simulation techniques, with better and more rigorous ways to approximate the quantum mechanical behavior of atoms and molecules. This improvement in connection between abinitio and MD calculations results has led to more accurate classical molecular dynamics (MD) simulations. 9 However, even atomistic analysis approach has its own length and time scale limitations. Due to explicit atomistic modeling methods, time scales are on the order of 1 fs and length scales on the order of 1 Å. Because of these very small time and length scales, typical atomistic simulations are limited to very small systems over very short times. Furthermore, there exist several multiscale methods to combine atomistic models and continuum models in a single simulation framework. These atomistic and multiscale methods have been successfully applied to investigate diverse defect structures within static or quasistatic descriptions. Moreover, even though one can model the material system with full atomic scale details, there are still some problems faced at mesoscale. As the number of atoms in a system increases, the possible local minimum energy configurations also grow rapidly (for N atom cluster the number of local minimum energy configurations grows faster than eN). Without knowing all the relative energy values of these local configurations, it is very difficult to prepare initial atomic configurations, which are most relevant to real physical systems. Secondly, the inaccuracy of an interatomic potential can introduce errors in relative energies of diverse configurations. Both of these problems seriously influence the reliability of atomistic simulations. Mesoplasticity is therefore needed to bridge this gap in between atomistic mechanisms of deformation and fracture to continuum behavior (refer Hartley 2003). Substantial progress has been made in the mesoscopic simulation area. The number of 3D dislocation loops required to represent fullscale plasticity of 10 submicron crystal is now manageable, and the solution requires integration of a few thousand equations of motion. However, the more challenging problem of polycrystalline plasticity will require additional breakthroughs, because of the conceptual difficulties of connecting dislocation dynamics (DD) to crystal plasticity models in a selfconsistent fashion. The main task is to retain the characteristic length scale from the discrete dislocation to the continuum description. Varieties of strain gradient theories have been proposed during the last decade. Most of them are formulated in a phenomenological manner and so the characteristic length is invoked into the theory, without rigorous treatment of its dislocation origins. In addition, it is also possible that the approach using a spectrum of length scales will be more appropriate to clearly establish the links between mesoscopic simulations and continuum mechanics. Recently, Ghoniem et al. (2002, 2003), Liu et al. (2004), Park and Liu (2004), and Lu and Kaxiras (2004) examined the efforts in multiscale simulation, that treat the problems involving either multiple scales in spatial or temporal or both simultaneously. Following sections give a brief review on the development of multiscale modeling and simulation of material behavior. In this review, we will be concerned with bridging length scales; while bridging time scales from the picoseconds times of atomistic vibrations to the micro, milli and larger time’s scales of various defect motions is an equally difficult but fundamentally different 11 matter. Moreover, this review mainly focuses on atomic atomistic/continuum coupling. Finally, some aspects along with recent work in coupled atomistic/continuum simulation through mesoplasticity will be reviewed in sections 2.3 and 2.4. 2.2 Coupled atomistic and continuum simulation approach Recently an alternative approach was explored that couples the atomistic and continuum approaches. The main idea behind this approach is to use atomistic modeling where the displacement field varies on an atomic scale and to use the continuum approach elsewhere. The critical part in these coupled simulations is the communication between the two descriptions of the materials. At the atomistic level, MD is the preferred mode while at the continuum level, FEM has been the dominant numerical technique used. Seamless coupling is required for the combined MD/FEM simulations, which is a key issue in multiscale simulations. The main difficulty lies in proper matching between the atomistic and the continuum regions. Since the details of the lattice vibrations, the phonons, which are an intrinsic part of the atomistic model, cannot be represented at the continuum level, conditions must be such that the phonons are not reflected at the atomisticcontinuum interface. Since the atomistic region is expected to be a very small part of the computational domain, violation of this condition quickly leads to local heating of the atomistic region and destroys the simulation. In addition, the matching at the atomisticcontinuum interface must be 12 such that largescale information is accurately transmitted in both directions. Some major multiscale atomistic/continuum simulations techniques developed are outlined below. 2.2.1 Finite element (FE) and molecular dynamics (MD) coupled simulation An early suggestion for coupling the continuum and the atomistic regimes was due to Kohlhoff et al. (1991), Hoover et al. (1992) and recently to RafiiTabar et al. (1998) and Abraham et al. (1998). In atomistic modeling of defects with long range stress fields, the boundary conditions applied to the atomistic region should be far from the defects. This results into unnecessary computations in the atomistic region far away from the defects. The other possible approach recognized was to surround the atomistic region with continuum region, which will provide flexible boundary conditions to the atomistic region. Seamless coupling is required for the combined MD/FE simulations, which is a key issue in multiscale simulations. Two methods used to combine the MD and FE region exist. One is through the imaginary surface, and the other is through the overlapping belt. The hybrid MD/FE method introduces a MD/FE imaginary surface as the handshake region to interface between the MD and the FE regions. MD atoms and FEnodes are coupled by mutual displacement boundary conditions. The kinetic and strain (potential) energy are defined for the entire system, including the handshaking interactions at the interface. The FE mesh in the MD/FE 13 handshaking region is scaled down to coincide with the atomic lattice. A onetoone correspondence of MD atoms and FEnodes is required in the interface region as shown in Fig. 2. Figure 2. Representative illustration of FE/MD handshake at the interface (Rudd and Broughton, 2000). In Fig. 2, The FE/MD interface is represented by a dashed vertical line. FE elements contributing to the overall Hamiltonian with full weight have dark shading, while those contributing with half weight have light shading. Two and threebody atomic interactions crossing the interface also carry half weight, and are shown with dotted lines. The FE mesh is gradually made coarser in the area FE region MD region MD/FE region MD/FE interface 14 far away from the interested atomistic region. Additionally, the width of the transition region is taken to be equal to the cutoff distance of the interaction potential used in the MD region. This provides a complete set of neighbors within the interaction range for all atoms in the MD region. Atoms/FEnodes that belong to the interface region not only interact via the interaction potential with the MD region but also are part of the nodes in the FE region. The positions and velocities of these atoms/FEnodes in the interface region (both on the MD side and on the FE side) must be consistent with each other. Equilibrium is guaranteed if the elastic constants in the continuum region match those of the atomistically modeled region. Figure 3 The stress waves propagating through the slab using a finely tuned potential energy color scale (Abraham et al. 2000) Kohlhoff et al. (1991) and Gumbsch et al. [1995 (a)] used the above approach in their ‘FEAt’ model for atomic scale fracture processes. They MD FE FE MD 15 obtained reasonably accurate calculation of critical energy release rates for fracture and crack tip blunting, allowing for a systematic comparison between the atomistic simulations and continuumbased criteria predicting brittle vs. ductile response. Fig. 3 shows results of coupled FE/MD simulation on rapid brittle fracture of a silicon slab flawed by a microcrack at its center and under uniaxial tension (Abraham et al. 2000). The stress wave propagation is seen in the slab through FE/MD interface with no visible reflection. Lidorikis et al. (2001) demonstrated a detailed comparison between hybrid simulation results and full MD multimillionatom simulations in the study on stress distributions in silicon/siliconnitride nanopixels. The coupled results showed a good agreement with full multimillionatom MD simulations. The atomic structures at the latticemismatched interface between amorphous silicon nitride and silicon induced inhomogeneous stress patterns in the substrate that cannot be reproduced by a continuum approach alone. However, the approach by Kohlhoff et al. (1991) does not support continuum dislocations in the finite element region, nor can it easily adapt to the size of the atomistic region during the simulation. The FEAt model requires that the atomistic region be specified at the start of the computation, thereby limiting any nonlinear deformation (including dislocations) to a predefined region The other method, which uses an overlapping belt, was found to be more flexible. In the overlapping area each node is surrounded by a group of atoms. 16 The forces exerted on the atoms in that region, due to the interaction with the MD region, make up the external forces on the FEnodes for the FE simulation. Tan and Yang (1994) and Yang et al. (1994) developed a combined continuum and atomistic simulation approach to simulate the transmission of crack tip dislocations from the atomistic assembly to the overlapping continuum. Noguchi and Furuya (1997) proposed a method that combines MD at the crack tip with a linear elastic outer domain in simulating elasticplastic crack advance. MD region was successfully moved as the crack propagated (Furuya and Noguchi, 1998, 2001). 2.2.2 Quasicontinuum method. Quasicontinuum (QC) theory developed by Tadmor et al. (1996a, 1996b) has been applied to a wide range of problems including dislocation structures, interactions of cracks with grain boundaries (Shenoy et al. 1998), nanoindentations (Shenoy et al. 2000, Smith et al. 2000, Smith et al. 2001), dislocation junctions (Rodney and Phillips 1999), atomistic scale fracture process (Miller et al. 1998). The results obtained are quite encouraging. Ortiz and Phillips (1999) reviewed the QC theory. Ortiz et al. (2001) established the same for multiscale simulations. QC seamlessly couples atomistic and continuum regions. In this approach the atomistic and finite element calculations are performed concurrently rather than in sequence, and therefore it is categorized as a concurrent multiscale approach. The QC theory initializes from a fundamental conventional atomistic model and systematically eliminates redundant degrees of 17 freedom by introducing kinematic constrains. The simulation initializes from atomistic scale and coarsens to describe the continuum region. While doing this, kinematic constraints are thoughtfully introduced by preserving the full atomistic description wherever required. Figure 4. Transition region for QC method (Curtin and Miller 2003) The QC method defines two types of atoms, ‘local representative atoms’ and ‘nonlocal representative atoms’, rather than identifying atomistic and continuum regions. In practice, however, the regions containing nonlocal representative atoms are essentially equivalent to the fully atomistic regions of other methods. Similarly, a local representative atom is coincident with either a continuum FE node or an atomic position near a Gauss point used to define the 18 energy of the continuum element. The language of nonlocal and local clearly associates the nonlocal atoms with ‘real’ atoms and the local atoms with the local FE region. The transition region between nonlocal and local representative atoms in the QC model is depicted in Fig. 4. The interface atoms are included in the atomistic energy but are also FE nodes of the continuum. The pad atom positions are dictated by interpolation from the FE nodal positions. The pad atom energies are not included but the energies of the real and interface atoms include their interactions with the pad atoms. To avoid overcounting of the energy of the interface atoms/nodes, the energies of the continuum elements adjacent to the interface (i.e. elements that have nodes corresponding to the interface atoms as shown in grey in Fig. 4) are weighted differently in the total potential energy sum. The total potential energy of the QC model is then obtained by summing the energies of all atoms in the atomistic region and at the interface and all elements in the continuum domain In this method, interatomic interactions are used in the model via the discrete calculation based on a local state of deformation. It starts by computing the energy of solid as a function of atomic positions. The configurations space of the solid is reduced to a subset of representative atoms. The over all equilibrium equations are then obtained by minimizing the potential energy described for the reduced configuration space. The total energy, in principle, needs to be calculated by summing the energy over all the atoms in system. N i 1 tot i E E , (2.1) 19 where N is the total number of atoms in the body. This full summation is bypassed by introducing approximate summation rules. The lattice quadrature analog of Eqn. (2.1) is written as Nr i 1 tot i i E n E where, ni is the quadrature weight that signifies how many atoms a given representative atom stands for in the description of the total energy, and Ei is the energy of ith representative atom which is summed over Nr representative atoms to get the total energy. In QC approach, displacement fields are obtained using FE method. The positions of the representative atoms are also calculated. From the new displacement field, again, energies are determined using atomistic calculations and the simulation continues. The method has been applied to single crystal FCC metals and shown to reproduce lattice statics results for a variety of line and surface defects (Tadmor et al. 1996b). Shenoy at al. (1998,1999) generalized the QC method and extended it to treat polycrystalline materials. As the first test of the model, they computed the grain boundary (GB) energy and atomic structures for various symmetric tilt GB’s in Au, Al, and Cu using both direct atomistic calculations and the model calculations. They found excellent agreement between the two sets of calculations, indicating the reliability of the model for their purpose. The generalized method was applied to a number of problems including the interaction of dislocations with grain boundaries in Al and the mechanics of steps on grain boundaries (Shenoy et al. 1998). In the study of Al, they used 20 nanoindentation induced dislocations to probe the interaction between dislocations and GB’s. Figure 5(a) FE mesh used to model dislocationGB interaction using QC method (Shenoy et al. 1999) Specifically, they considered a block oriented so that the (111) planes are positioned to allow for the emergence of dislocations which then travel to the GB located at ~200Å beneath the surface. Fig 5 (a) shows FE mesh for nanoindentation model used for multiscale simulation using QC method. The surface marked AB is rigidly indented in order to generate dislocations at A. First, the energy minimization is performed to obtain the equilibrium configuration of the GB, and then the FE mesh is constructed accordingly. Full atomistic description is used for the region that is expected to participate in the dislocation 21 GB interaction, while in the far fields the mesh is coarser. The displacement boundary conditions at the indentation surface are then imposed onto this model structure, and after the critical displacement level is reached, dislocations are nucleated at the surface. They use the EAM potential to calculate the energies at atomistic level. They observed closely spaced (15 Å) Shockley partials nucleated at the free surface. Figure 5(b) Snapshots of atomic positions at different stages in the deformation history. (Shenoy et al. 1999) As seen in [Fig 5(b)], the partials are subsequently absorbed at the GB with the creation of a step at the GB and no evidence of slip transmission into the adjacent grain is observed. The resultant structure can be understood based on the concept of displacement shift complete lattice associated with this symmetric tilt GB. As the load is increased, the second pair of Shockley partials is 22 nucleated. These partials are not immediately absorbed into the GB, but instead form a pileup [Fig. 5(b)]. The dislocations are not absorbed until a much higher load level is reached. Even after the second set of partial dislocations is absorbed at the GB, there is no evidence of slip transmission into the adjacent grain, although the GB becomes much less ordered. The authors argued that their results give hints for the general mechanism governing the absorption and transmission of dislocations at GB’s. In the same work they also studied the interaction between a brittle crack and a GB and observed stressinduced GB motion (the combination of GB sliding and migration). In this example, the reduction in the computational effort associated with the quasicontinuum thinning of degrees of freedom is significant. For example, the number of degrees of freedom associated with the mesh is about 104, which is three orders of magnitude smaller than what would be required by a full atomistic simulation (107 degrees of freedom). Miller et al. (1998) extended the method to study fracture mechanics at the atomic scale. Fig. 6 shows a crack tip subjected to mode I loading approaching a grain boundary in a nickel bicrystal, with increasing levels of external load. The colors indicate levels of slip from zero (red) to maximum (blue). The snapshots show initial dislocation from the grain boundary (t=3), followed by crack extension by cleavage (t=4) and in the end grain boundary migration to the crack tip (t=6). They also investigated the effect of grain orientation on fracture and the interaction of cracks with grain boundaries. 23 Rodney and Phillips (1999) extended the method to three dimensions and used it to study junction formation and destruction between interacting dislocations. Further, the method was also extended to treat complex Bravais lattice materials (Tadmor et al. 1999). Figure 6 Quasicontinuum simulation of fracture at atomic scale (Miller et al. 1998) The authors refer to using the “local” regime of the QC method, which refers to the case where each FE node represents a very large number of 24 atomistic degrees of freedom, which is modeled as corresponding to an infinite solid homogeneously deformed according to the local strain at the node. In this limit, it is advantageous to employ effective Hamiltonians to compute energetics for the representative atoms. Such Hamiltonians can be constructed from abinitio calculations, and may include physics that atomistic simulations based on classical interatomic potentials (such as EAM) are not able to capture. For example, by constructing, an effective Hamiltonian parameterized from ab initio calculations, Tadmor and his group was able to study polarization switching in piezoelectric material PbTiO3 in the context of the quasicontinuum model in the local limit (Tadmor et al. 2002). With this effective Hamiltonian, it was shown that the behavior of a largestrain ferroelectric actuator under the application of competing external stress and electric fields could be modeled successfully, reproducing, for example, all the important features of the experimental strain vs. electric field curve for the actuator. The advantage of these simulations is that they also provide insight into the microscopic mechanisms responsible for the macroscopic behavior, making possible the improvement of design of the technologically useful materials. QC was also used to study nanoindentation into silicon single crystals (Smith et al. 2000, 2001) and of aluminum thin film (Tadmor et al. 1999). Fig. 7 shows nanoindentation simulation of a thin film aluminum that was carried out. Dislocations are seen emitted from indenter corners (indenter not in the picture). 25 Figure 7 Nanoindentation simulation on thin film Al using QC method (Tadmor et al. 1999) Knap and Ortiz (2001) presented a streamlined and fully nonlocal threedimensional version of the QCtheory. However in spit of all this, the QC method does not allow for the continuum description of a dislocation and therefore every dislocation in the QC model requires a full complement of atomistic degrees of freedom around its core and slip plane. This means that only a few dislocations can be present in a QC simulation before the computational cost becomes comparable to that of the fully atomistic description. Additionally the computational effort in the QC approach quickly approaches that of the fully atomistic model once a relatively small number (less than 100) of dislocations are developed in the microstructure. 26 2.2.3 Lattice Green’s function. Thomson et al. (1992), Zhou et al. (1993) and SchiØtz et al. (1997) used this multiscale approach for calculating the lattice defects. In this approach, the actual number of atomic degrees of freedom is not changed. They treat a large portion of atomic response as linear instead of nonlinear while they are still explicitly modeled. As a result they considerably reduced the computational time required. Rao et al. (1999) extended the Green’s function technique to two and three dimensions. They demonstrated this technique by applying it to relax the boundary forces in the simulations of crossslipped core structures of (a/2)[110] screw dislocations of FCC structure. Further, MasudaJindo et al. (2001) used this approach to study the mechanical properties of nanoscale materials. The initial atomic structures of the dislocations and cracks were determined both from the elastic solutions as well as from those by lattice Green’s function method for the infinite systems. They calculated the Green function for the defective lattice, with dislocation and crack, by solving the Dyson equation, appropriate for absolute zero temperature. The thermal expansion and the temperature dependence of the interatomic force constants are determined using the statistical moment method and they are taken into account in the lattice Green’s functions. With this methodology, they investigated the strength and fracture properties for the nanocrystalline materials like semiconductor quantum wire and carbon related materials like graphenes and nanotubes. The O(N) tightbinding molecular dynamics (TBMD) method is 27 used to analyze the reconstruction of atomic bonding near the crack tip as well as the cleaved surface. . Overall, Lattice Green’s function approach is used for static systems and is effective in seeking longterm balanced configurations. 2.2.4 Coarsegrained molecular dynamics. Rudd and Broughton (1998) took into account the limitation that the continuum based origins of the FE method lead to difficulties in a smooth transition from the atomistic to the continuum. Thus, they reformulated the continuum region using a method called coarsegrain molecular dynamics (CGMD), in which they developed the degreeoffreedom reduction as a more natural extension of the underlying discrete molecular dynamics. CGMD substitutes FE mesh used outside the atomistic region, and it connects to MD at the atomistic scale. In this approach, the equations of motion do not rely on the assumptions of continuum elasticity theory. As a result, CGMD was found more appropriate for mesoscopic elastic solids. The model aims at constructing scaledependent constitutive equations for different regions in a material. In their methodology, the material of interest is partitioned into cells, whose size varies so that in important regions a mesh node is assigned to each equilibrium atomic position, whereas in other regions the cells contain many atoms and the nodes need not coincide with atomic sites. The CGMD approach produces equations of motion for a mean displacement field defined at the nodes by defining a conserved energy functional for the coarsegrained system as a 28 constrained ensemble average of the atomistic energy under fixed thermodynamic conditions. The key point of this effective model is that the equations of motion for the nodal (mean) fields are not derived from the continuum model, but from the underlying atomistic model. The nodal fields represent the average properties of the corresponding atoms, and equations of motion (in this particular case Hamilton’s equations) are constructed to describe the mean behavior of underlying atoms that have been integrated out. One important principle of CGMD is that the classical ensemble must obey the constraint that the position and momenta of atoms are consistent with the mean displacement and momentum fields. When the mesh nodes and the atomic sites are identical, the CGMD equations of motion agree with the atomistic equations of motion. As the mesh size increases some shortwavelength degrees of freedom are not supported by the coarse mesh. However, these degrees of freedom are not neglected entirely, because their thermodynamic average effect has been retained. This approximation is expected to be good if the system is initially in thermal equilibrium, and the missing degrees of freedom only produce adiabatic changes to the system. The Hamiltonian was derived originally for a monoatomic harmonic solid, but can be easily generalized to polyatomic solids (Rudd and Broughton 1998). After deriving the equations of motion from the assumed Hamiltonian for a particular system, one can perform the CGMD for the nodal points. As a proof of principle, the method was applied to onedimensional chains of atoms with periodic boundary conditions where it was shown that the CGMD 29 gives better results for the phonon spectrum of the model system compared to two different FE methods (Rudd and Broughton 1998). A variety of other calculations have also been performed with the CGMD to validate its effectiveness (Rudd and Broughton 2000, unpublished). Rudd (2001) developed a nonconservative version of CGMD to remove unphysical elastic wave reflection. Further, Rudd (2002) used CGMD to study the behavior of submicron MicroElectroMechanical Systems (MEMS), especially microresonators, often called NanoElectroMechanical Systems (NEMS). Although the CGMD has proven to be reliable in the description of lattice statics and dynamics, the implementation of the model varies from system to system. This is because different approximations have to be made to the Hamiltonian that represents a particular system. On the other hand, such approximations can be estimated and controlled in the CGMD method. This advantage makes the CGMD method a good substitute for the FE method in the MAAD approach when a high quality of FE/MD coupling is required. The CGMD approach resembles the quasicontinuum model in the sense that both approaches adopt an effective field model, an idea rooted in the renormalization group theory for critical phenomena. In both approaches, less important (long wavelength) degrees of freedom are removed while the effective Hamiltonian is derived from the underlying finescale (atomistic) model. Although both approaches are developed to couple FE and atomistic models, the quasicontinuum method is mainly applicable to zero temperature calculations while the CGMD is designed for finite temperature dynamics. The 30 quasicontinuum model has shown its success in many applications, but the CGMD approach has yet to show its wider applicability and versatility. The CGMD results obtained were limited to validation problems that demonstrated the capabilities of the method. 2.2.5 Compound wavelet matrix method. This method was found suitable for multiscale investigations. Frantziskonis and Deymier (2000) designed compound wavelet matrix (CWM) method to bridge multiscale models over different ranges of spatial and temporal scales. They applied CWM to analyze the microstructure of 2D polycrystalline systems as well as their evolution during normal grain growth. On the finer scale MD simulation with LennardJones (LJ) system, and on the coarser scale MC simulation of the Qstates Potts model was used. The basic idea used was to produce matrices of the wavelet coefficient from energy maps representing the spatial distribution of the local excess energy in the microstructures obtained with both methodologies. The full description of the material was then obtained by merging the matrices of the wavelet coefficients representing the material at different scales through CWM method. The CWM then characterizes the materials over a range of scales that is the union of the scales treated by the two methodologies. This method possesses several advantages. First, it does not assume a priori that a collection of small microscale systems is equivalent to a microscalebased model of a large system. 31 Second, the simulation time of the coarsest methodology is not controlled by the methodology with the slowest dynamics. Figure 8 Schematic diagram of the construction of CWM (a) from the LJ system of 1x1 units of length discretized into 256x256 points (b) from the Potts system of 2x2 units of length discretized into 128x128 points. (Frantziskonis and Deymier, 2000) To present an illustrative example of the CWM method applied to 2D grain growth, an MD simulation of a 2D LJ system and a MC simulation of a Qstates Potts model that can overlap over a range of spatial and time scales were presented, Fig. 8. The bridging of these simulations was accomplished by the 32 overlap in scale of the ‘mesoscopic’ Qstates Potts model with the atomistic LJ model and constructing the CWM. In Fig. 8 the arrows indicate the substitution (transfer of matrix statistics) of matrices in forming the compound matrix. Only three such substitutions are shown for illustration (Frantziskonis and Deymier, 2000). 2.2.6 Other approaches. Another recent development in coupled method includes algorithms to pass defects between that atomistic and continuum regions developed by Shilkrot et al. 2002 and given the name, the coupled atomistics and discrete dislocation (CADD) method. The CADD approach is similar to the FEAt approach, in that it begins from a separation of the physical problem into spatial regions that are modeled either by a fully atomistic description or by continuum finite elements. CADD extends these methods, however, by allowing for the presence and movement of continuum discrete dislocations in the continuum regime that interact with each other and with the atoms in the atomistic region via their elastic fields. The philosophy behind the CADD method is that atomistic resolution is required in some regions of a material, but away from those regions the deformation is well described by elasticity plus the plasticity due to the motion of welldefined continuum dislocations whose interactions can be described by elastic continuum fields. 33 Some other concurrent approaches are similar to the MAAD method but concentrate on linking two different length scales rather than three. For example, Bernstein and Hess (2001) have simulated brittle fracture of Si by dynamically coupling empiricalpotential MD and semi classical TB approaches. 2.2.7 Lattice material point method (LMPM) 2.2.7.1 Material Point Method (MPM). A new computational method, namely, the material point method (MPM) was developed by Sulsky and others (1995,1996) at the University of New Mexico with support from the Sandia National Laboratories. MPM was developed from FLIP (Fluid  implicit  particle) particle – in – cell (PIC) method by extending the computational fluid dynamics capability of the code, to solving solid mechanics problems. Sulsky et al. (1994), initially considered application of MPM to 2 dimensional impact problems and demonstrated the potential of MPM. Sulsky and Schreyer (1996) gave a more general description of MPM, along with special considerations relevant to axisymmetric problems. The method utilizes a material or Lagrangian mesh defined on the body under investigation, and a spatial or Eulerian mesh defined over the computational domain. The set of material points making up the material mesh is tracked throughout the deformation history of the body and these points carry with them a representation of the solution in a 34 Lagrangian frame. Interactions among these material points are computed by projecting information they carry onto a background finite element mesh where equations of motion are solved. They reported that the material point method does not exhibit locking or an overly stiff response in simulations of upsetting. In MPM, a solid body is discretized into a collection of material points, which together represent the body. As the dynamic analysis proceeds, the solution is tracked on the material points by updating all required properties, such as position, velocity, acceleration, stress, and temperature. At each time step, the material point information is extrapolated to a background grid, which serves as a computational scratch pad to solve the equations of motions. The solution on the gridnodes is transferred using finite element shape functions, on the material points, to update their values. This combination of Lagrangian and Eulerian method has proven useful for solving solid mechanics problems involving materials with history dependent properties such as plasticity or viscoelasticity effects. MPM is amendable to parallel computation (Parker 2002), implicit integration method for unconditionally stable time increments (Guilkey and Weiss 2003, Sulsky and Kaul 2003, Cummins and Brackbill 2001) and alternative interpolation schemes using C1 continuous interpolating functions for improved smoothness, without numerical noise, in field properties (Bardenhagen and Kobar 2000). The computational methodologies of material point method as well as of molecular dynamics are discussed in Chapter 4. 35 2.2.7.2 Coupled MPM/MD coupling approach Tan (2001) proposed this method. The basis for LMPM, coupling MD and MPM, is that a continuum material point is also an aggregate of atoms. MD/MPM connection was realized by the Lattice Material Point (LMP) and the background grid. LMP has two viewpoints  a lattice from an atomistic perspective and a material point from a continuum perspective. Providing a smooth transition or handshake between these two representations is the spirit of LMPM. From an atomistic point of view, a LMP represents a lattice, a set of atoms. Each LMP has corner points to define the occupied space. The corner point updates the position according to the interpolated velocity field from the MPM grid. The LMP continuum stresses are the average of atomic stresses, which are evaluated through statistical average of the force on each atom from the local neighboring atoms. LMPM uses a uniform background grid that hosts both continuum cells and atomistic cells. The simulation region is divided into two regions  an atomistic region and a continuum region; grid nodes connect the two regions. Nodes at the interface between the continuum region and atomistic region, such as node m in Fig. 9, collect stress and temperature information from LMP in both continuum and atomic sides. On the continuum side, the continuum LMP (gray points) interpolates the value from grid nodes and updates the values at material points (black points). At atomic side, the MD simulation needs the neighboring atomic position information to compute the potential and the forces on each atom (black points). 36 Figure 9. Connection between continuum region and atomistic region (Tan 2001). Using LMPM Tan (2001) simulated an intergranular crack that propagates to an intergranular SiC under mode I extensional load. LMP was used to convert atomistic dislocations to continuum ones. When an atomistic dislocation moves to the continuum region, it was substituted in situ by a dislocation of the same Burgers vector but embedded in the MPM described continuum. As the consequence of the newly created continuum dislocation, the extra halfarray of atoms in the atomistic assembly were removed from the MD calculation to retain the global conservation of mass. In the continuum region, dislocation motion was directed by the crystal slip system with speed given by the dislocation dynamics curve. 37 Ayton et al. (2001) developed a computational methodology, which included the continuum level simulation using MPM and the microscopic method of nonequilibrium MD. They examined the behavior of a membrane bilayer/solvent system using this microtomacro dynamical feedback simulation technique. With this method, two simulations at different time and length scales were interfaced into a unified simulation methodology. The interface was accomplished via an information transfer where selected material properties (transport coefficients) and state parameters (density) were calculated in one spatial/temporal regime and then used as initial input in another leading to a closed feedback loop. The multiscale methods based on Lattice Green’s function and QC methods are static methods. These give the longterm balanced configurations. Also time scale is not an issue in these approaches. However dynamic effects are important at nanoscale, (refer Tan 2003). Dynamic analysis at nanoscale is possible using coupled FE/MD and LMPM methods. The material point/atoms linkage in LMPM has some similarity with CGMD, where a material point can be looked as a coarsegrained atom group. However, in FE/MD coupling, one—toone correspondence of FEnode and MDatom provides extra constraint and blocks the moving of atomic defects into continuum region (Tan 2001). Even in CGMD this may become a difficult issue. Additionally, under large deformations, the FE elements become distorted and it is difficult to implement the entire computational process in the FE simulations. 38 2.3 Multiscale modeling of plasticity using coupled mesoplasticity and continuum simulation. Two complimentary approaches in the area of mesoscale analysis have been advanced. First of the two approaches involves dislocation dynamics and the second approach is using continuously dislocated continuum theory or statistical mechanics approach. Recently numerous efforts have been devoted for coupled mesoscopic and continuum scale simulations. Herein the work contributed by Zbib and Rubia (2002), Hartley (2001, 2003), Khan et al. (2004),, Abu AlRub and Voyiadjis (2004), LeSar and his coworkers is reviewed. Dislocation Dynamics (DD) has been proven useful for simulating material behavior at mesoscale. In this strategy, instead of simulating the dynamics of atomic systems, the dynamics of defect ensembles are simulated in the material. Examples of this approach are the dynamic simulations of interacting cracks in brittle materials, or dislocations in crystalline materials. DD theory involves incorporating nonlinear interaction among dislocations and interaction of dislocations with interfaces. Lepinoux and Kubin (1987) were the first to develop twodimensional models of dislocations, although an alternate version. More developed versions were then proposed (Ghoniem and Amodeo 1988a,b, Amodeo and Ghoniem 1990a,b, Gulluoglu et al. 1989, Lubarda et al. 1993, Barts and Carlsson 1995, Wang and LeSar 1995). During this time the modeling 39 capabilities were limited to twodimensional and consisted of infinitely long and straight models of dislocation in an isotropic infinite elastic medium. Using DD, complex phenomena such as dislocation cells, slip bands, microshear bands, persistent slip bands and dislocation tangles, were proved to dictate the material properties. To understand more realistic features of the microstructure in crystalline solids, Kubin and coworkers (Kubin and Canova 1992, DeVincre and Kubin 1994), pioneered the development of 3D lattice dislocation dynamics. More recent advances in this area have contributed to its rapid development (Zbib et al. 1998, Ghoniem et al. 2000, 2001). Another complimentary approach for mesoscale analysis has been based on a continuously dislocated continuum theory or the statistical mechanics approach (refer Needleman 2001) and is also termed as mesoplasticity. In this approach, models involving interactions of dislocation ensembles and the manner they influence material properties are included in the evolution equations based on classical continuum mechanics framework using representative volume element (RVE). These governing equations along with boundary conditions are solved for the complete description of the deformation problem. The mesoplastic constitutive relation was first proposed by Taylor (1938). Kröner (1970) proposed a theory of global plastic multiaxial deformation that incorporated a description of the local dislocation structure in the form of state variables. Mesoplastic constitutive relation has been refined and improved by many researchers, involving Hill and Rice (1972), and Asaro and Needleman (1985). With the 40 increasing computational power, the relation can be implemented in some FEM codes to solve complex threedimensional problems. For example, using the mesoplastic theory, Peirce et al. (1982) analyzed numerically the nonuniform and localized deformations of ductile single crystals subjected to tensile loading. Huang (1991) developed a usermaterial subroutine incorporating mesoplasticity in the commercially available ABAQUS FE implicit program. This twodimensional subroutine has implemented the theoretical framework of Hill and Rice (1972). Kalidindi et al. (1992) developed an implicit timeintegration procedure mainly based on the constitutive model of Asaro and Needleman (1985). Kalidindi et al. simulated the loaddisplacement curve in a planestrain forging on copper that captured successfully major features of the experimental data. They reported jumps on the calculated curve although no jumps were observed from the experimental data. Yoshino et al. (2001) applied the mesoplasticity theory in 2D FEM to simulate the dislocation generation and propagation during indentation of single crystal silicon. Kysar (2001) studied the crack growth along a copper/sapphire bicrystal interface. Wang et al. (2004) studied the dependence of nanoindentation pileup patterns and microtextures on the crystallographic orientation on single crystal copper using a conical indenter. They, however, reported an order of magnitude’s deviation of the numerical loaddisplacement results from experimental data. Fivel et al. (1998) developed a 3D model to combine discrete dislocation with FEM for nanoindentation simulation. 41 Above mentioned models are based on basic deformation mechanisms and aim at predicting strength properties of metals at small length scales. However sizedependent problems which include, modeling of dislocation boundaries, dislocation in heterogeneous structures like multiple layer and thin film structures, dislocation interaction with interfaces and associated shape changes with lattice rotations, and additionally deformation in nanostructured materials, localized deformation and shear bands still cannot be addressed using the technique of DD alone due to computational limitations. Hence coupled DD and continuum scale approach is being developed. The main difficulty in these models, among other things, lies in connecting the observed parameters (used in continuum models), such as stress and strain, to the collective behavior of dislocation groups and their interaction with particles and interfaces. While within the continuum mechanics framework, the governing equations of the material response are developed based on a representative volume element (RVE) over which the deformation field is assumed homogeneous, the DD models are presumably capable to describing the heterogeneous nature of the deformation field within the RVE. Thus with a proper homogenization theory, or field averaging, one can couple the discrete DD models with continuum approaches and provide a rigorous framework for analyzing deformations at small scales, where surfaces and interfaces are of great importance in determining material response. van der Giessen and Needleman (1995) developed this approach. They coupled a twodimensional discrete dislocation dynamic model with continuum finite element model. Further, 42 Zbib and Rubia (2002) extended the framework to merge the two scales, namely nanomicroscale, where plasticity is determined by explicit threedimensional dislocation dynamics (DD) analysis, and the continuum scale, where energy transport is based on basic continuum mechanics laws. The output of this approach is a coupled threedimensional discrete DD with continuum elastoviscoplastic model. A method based on superposition principle is adopted for modeling dislocation surface interactions in both homogeneous and heterogeneous materials. Finite element method (FEM) is the technique adopted for continuum scale analysis. Three models are combined together, the DD model, the solid mechanics and the heat transfer model. Their hybrid multiscale model of plasticity is capable of explicitly modeling complex smallscale plasticity phenomena, including surface effects, ledges, microshear bands, dislocation boundaries, deformation of thin layers. Their results on multiscale nanoindentation simulation, for investigating deformation and strength in thinlayered structures are shown in Fig. 10 and 11. The test sample used was Fe3% Si single crystal and Oxide thin film of thickness 10.0 nm on it and the load was applied using a Berkovich indenter (192 nm tip radius of curvature) and in a direction normal to (001) plane. For more details, refer to Zbib and Rubia (2002). 43 Figure 10 Nanoindentation simulation stress contour using multiscale model of plasticity. (Zbib and Rubia 2002) The stress field using multiscale nanoindentation simulation is shown Fig 10. For the case of random distribution of dislocation sources the evolution of dislocations is shown in Fig 11 (a) and the depth during nanoindentation as predicted by their hybrid multiscale model is shown in Fig 11(b). The authors comment that the resulting dislocation structure and predicted depth are consistent with the experimental observations reported by Zielinski et al. (1995) and Bahr and Zbib (unpublished). 44 Figure 11 Results of nanoindentation simulations using multiscale plasticity model. (a) developed dislocation structure underneath the indenter. (b) Predicted indentation depth showing the development of surface distortions. (Zbib and Rubia 2002). Further, Khan, Zbib and Hughes used multiscale DD plasticity to model planar dislocation boundaries in deformed FCC single crystals. For multiscale simulations, they constructed a pure tilt boundary and experimentally observed extended geometrically necessary boundaries (GNBs) within the representative volume element (RVE). Coupling between discrete DD analysis approach with 45 continuum finite element is obtained to correct for the boundary conditions and image stress. Their approach is that, they calculate the PeachKoehler force acting on a dislocation segment from the stress field of all other dislocations and the applied stresses. And this is calculated while the plastic deformation of the single crystal is obtained explicitly accounting for the evolution of a multitude of dislocation loops and curves. Then following the standard finite element procedure and using linear interpolation shape functions over the segment, the PeachKoehler force per unit length is distributed equally to the nodes. Thus, once all the forces are assembled, the net force on each node would have contributions from all the segments connected to it. The Dislocation nodes move simultaneously in the glide direction over a characteristic time corresponding to the least time increment required for an interaction to take place. Since the governing ewuation of glide motion for each dislocation node is nonlinear, the motion and interaction of an ensemble of dislocations in a 3D crystal is integrated over time yielding macroscopic plastic distortion. The authors conclude that the right boundary condition of the RVE is critical in modeling GNBs and their longrange stresses. They also considered the effects of various numerical factors such as domain length and mesh sensitivity. Additionally the effect of changing the spacing between two dislocation boundaries on the selfstress field and the stability, particularly in the space between the two dislocation boundaries, is also discussed. 46 Towards an effort to couple the classical continuum mechanics with mesoplasticity, Hartley (1985, 2001) proposed the concept of a dislocation mobility tensor. Recently Hartley (2003) developed a method for linking thermally activated dislocation mechanisms of yielding with continuously dislocated continuum plasticity theory as developed by Kröner (1958) and Bilby (1955). For information transfer in coupled simulations, he proposed that the total distortion needs to be decomposed into components that are due to dislocation motion alone and to the subsequent lattice distortion caused by the dislocations and external boundary conditions. The volumetric distortion associated with movement of dislocations should be used to connect to continuum theories. To attain the same, he introduced two quantities, the dislocation mobility tensor and the dislocation velocity vector, that characterizes the geometrical nature and stress dependence of dislocation ensembles in a way that permits their inclusions in continuum theories of elastoplastic deformation. The magnitude of the dislocation density vector for a slip system describes the net length, or volume average, of dislocation lines having the same Burgers vector in a RVE. The orientation of the vector is determined by the relative proportion of the total length due to screw and edge components. The vector is constructed by first defining a dislocation distribution vector for an arbitrary direction as the product of a unit vector in the direction and the total length per unit volume of dislocation segments parallel to that direction. Components of the dislocation density vector relative to a reference coordinate system are obtained by integrating the corresponding components of the dislocation distribution vector over all possible 47 orientations. The magnitude of the screw component and the orientation of the vector relative to the slip direction are determined by projecting the dislocation density vector parallel to the slip direction. Defining the dislocation density vector in this manner permits the definition of the net Peach–Koehler force per unit volume as the virtual force on the dislocation configuration due to internal and external sources of stress. It is important to note that the dislocation density vector defined in this fashion only describes the screw–edge character of the net dislocation population. It does not contain details of the local geometry of the dislocation arrangement, such as the relative positions of dislocation segments; hence, it cannot be used for calculations that depend on these values. Above concepts can be applied directly to the kinematics and kinetics of dislocation dynamics calculations that produce detailed pictures of the dislocation arrangement in a computational volume. The description of dislocation density can be employed to insert the results of these calculations into continuum descriptions of macroscopic material behaviour in order that the dislocation dynamics calculations can be employed to derive constitutive equations applicable to macroscopic structures. The dislocation mobility tensor provides a convenient form for the inclusion of specific dislocationbased microscopic constitutive relations into a model that can then be used as the basis for dislocation dynamics calculations or further analytical treatments. Constitutive equations that describe the relationship of the applied stress to the associated dislocation strain rate at the onset of observable permanent deformation can be developed in terms of internal structure of the material from a model based on 48 thermal activation of dislocations over barriers in their glide plane. Using this approach Hartley demonstrate how stress, temperature and material parameters enter into an expression for the flow law using the intersection of a gliding dislocation with a forest dislocation as an illustration. Recently, Abu AlRub and Voyiadjis (2004) made use of micro to nanoindentation experiments to determine analytically the material intrinsic length scale of strain gradient plasticity theory. The enhanced gradient plasticity theories formulate a constitutive framework on the continuum level that is used to bridge the gap between the micromechanical plasticity and the classical continuum plasticity. They are successful in explaining the size effects encountered in many micro and nanoadvanced technologies due to the incorporation of an intrinsic material length parameter into the constitutive modeling. However, the full utility of the gradient type theories hinges on one’s ability to determine the intrinsic material length that scales with strain gradients, and this study aims at addressing and remedying this situation. Based on the Taylor’s hardening law, a micromechanical model that assesses a nonlinear coupling between the statistically stored dislocations (SSDs) and geometrically necessary dislocations (GNDs) is used in order to derive an analytical form for the deformationgradientrelated intrinsic lengthscale parameter in terms of measurable microstructural physical parameters. They also present a method for identifying the lengthscale parameter from micro and nanoindentation experiments using both spherical and pyramidal indenters. The deviation of the 49 Nix and Gao (1998) and Swadener et al. (2002) indentation size effect (ISE) models’ predictions, from hardness results at small depths for the case of conical indenters and at small diameters for the case of spherical indenters, respectively, is largely corrected by incorporating an interaction coefficient. It thus compensates for the proper coupling between the SSDs and GNDs during indentation. They also present experimental results which show that the ISE for pyramidal and spherical indenters can be correlated successfully by using the proposed model. LeSar and coworkers (LeSar and Rickman 2003) have devoted some efforts in developing coarsegraining strategies due to sizelimitation imposed by direct DD simulations. They have begun to tackle this issue by investigations of a number of specific problems, namely the combined motion of dislocations and solutes, and the problem of the longrange interactions of blocks of dislocation ensembles. In the limit where the solute mobility is high relative to that of dislocations, dislocation–solute atmosphere pairs can be treated as quasiparticles, and the degrees of freedom associated with individual solute atoms can be eliminated (Rickman et al. 2003). Thus, an effective mobility can be derived and can be related to the details of such interaction. However, in the other extreme limit of low solute mobility compared with dislocations, the solute atoms function essentially as static traps for dislocations. In this case, LeSar and coworkers modeled slip as a discrete onedimensional ‘birthanddeath’ stochastic process. These examples show that it may be possible to average out fast time 50 variables in DD simulations. The issue of spatial averaging has been recently dealt with by LeSar and Rickman (2002) and Wang et al. (2003b). Starting from the work of Kossevich (1979) on the interaction energy of systems of dislocations, an energy expression in terms of the dislocation density tensor (Rickman et al. 2003) is derived. The basic idea in this approach is to divide space into small averaging volumes to calculate the local multipole moments of the dislocation microstructure and then to write the energy (LeSar and Rickman 2002) or the force (Wang et al. 2003b) as an expansion over the multipoles. Numerical implementation of the multipole expansion for the Peach–Kohler force on a test dislocation separated from a volume containing a dense dislocation aggregate revealed that such a technique is highly accurate and can result in significant coarse graining in current DD simulations. 2.4 Coupled Atomistic and Discrete Dislocation plasticity Plastic deformation and fracture of ductile materials involves important physical phenomena at multiple length scales. Some phenomena (e.g., dislocation nucleation, mobility, crossslip, crack formation, and growth) are intrinsically atomistic. Atomistic studies can address these unit processes involving a few defects but are usually unable to address largerscale deformation except with supercomputers. Plastic deformation at the micron scale, due to motion and interaction of many dislocations, can be described by representing the dislocation core structure, which is the main idea used behind discrete dislocation dynamic modeling of material. A true coupled multiscale 51 model should be capable of dealing with the necessary atomistic degrees of freedom, including all longrange interactions, within one overall computational scheme. Key and distinguishable features of such a model are that, dislocations exist in the continuum region and that they are mechanically coupled to one another and to the atomistic region. Also such dislocations can be passed between the atomistic and the continuum regions so that the plastic deformation is not confined to the atomistic region. Passing of dislocations requires two main developments: (1) detection of dislocations near the atomic/continuum interface that are selected for passing and (2) a procedure for moving such dislocations across the interface. Shilkrot et al. developed such a method, which they call it as: the coupled atomistic and discrete dislocation (CADD) method. In order to detect the dislocations and pass them to the continuum region, they introduce a “detection band” of triangular elements between atoms, with each element sitting on three different slip planes, as shown in Fig. 12. 52 Figure 12. (a) CADD atomistic/continuum interface (c) One detection element, indicating the 3 slip planes on which it lies (Curtin and Miller 2003) If the dislocation passes along one of these planes, it generates a Lagrangian finite strain in the element, element, which can be used to identify the Burgers vector and slip plane of the nucleated defect. For a given crystal structure and orientation, all the possible slip systems and associated strains are known. The detection algorithm monitors the strains in the detection band elements and 53 compares them to the known dislocation slip strains after each energy minimization step. When the strain in an element corresponds to a particular dislocation, the dislocation core is assumed to reside at the centroid of that element. Use of lagrangian finite strain tensor allows for large deformation and rotations. Because multiple dislocations can pass through one element, it is necessary to keep track of all previous slip activity and consider only the displacements due to new defects. After this recognition step the detection establishes the location of the core, the slip plane, and the Burgers vector of the continuum entity. Then the core is artificially shifted along its slip plane from its location in the detection band to a location across the interface in the continuum region by adding the displacement fields associated with a dislocation dipole. These displacements cancel the original core in the atomistic region and add a new core in the continuum region. Once this core is in the continuum region, it is added to the array of discrete dislocation with care being taken to define the branch cut in the continuum displacement field such that it matches the slip plane along which the dislocation moves. CADD is well suited to problems including nanoindentation, fracture, grain boundary sliding, slip transmission across grain boundaries and void nucleation and growth in crystals. So far, this approach has only been implemented in 2D systems, and it has been shown to agree well with the 2D atomistic calculations for Al. Figure 13(a) shows the sample geometry, the rigid indenter and the fully atomistic region as a small box located near the indenter corner, where dislocation nucleation will occur. 54 Figure 13. (a) Nanoindentation simulation using CADD (b) Predicted forcedisplacement curve (c) Positions of nucleated dislocations (Curtin and Miller 2003) The CADD problem contains only about 3000 atoms/nodes (2800 atoms and 200 FE nodes) whereas the full atomistic problem contains about 100 000 atoms. Under indentation, dislocations are nucleated, move away from the indenter into 55 the continuum, and exert a back stress on the atoms near the indenter corner. The predicted forcedisplacement behavior and the positions of nucleated dislocations are shown in figures 13(b) and (c), and compare favorably to the atomistic results and are also presented. The differences in dislocation position are largely associated with the very shallow energy minima for the dislocations in the continuum region, which leads to relatively large differences in position caused by small errors in the predicted stress fields. This approach was also applied to study growth of an atomic scale void in FCC Ni, due to plastic deformation. Recently there have been significant modifications in the nature of atomistic/continuum coupling at the interface (Shilkrot et al. 2004) due to original method showing spurious effects. Compared to CADD the QC method requires full atomic resolution wherever the CADD discrete dislocations are present. 2.5 Conclusion Currently there appears to be a great need to develop general mathematical and computational methods for a truly seamless multiscale approach for computer simulations of nano and Microsystems. In these early stages of development, the computational techniques are developed within a specified range of space and time scales. It is understood that the transition between one spacetime range to another is carried out by a process of handshaking, that is the information gained from a lower scale is summarized into a finite set of parameters, and passed on to the higher scale. This is 56 acceptable as long as the parameters represent a rigorous reduction of the enormous DOF of a lower length scale into a few generalized DOF. However, the theoretical foundations and computational implementation of a more rigorous process remains unresolved. The smooth transitions between various length scales can be ensured is there could be generalizations of the concept of DOF from those associated with spacetime (e.g. geometry), to those representative of statistical configurations (e.g. conductivity, mobility, etc.). 57 CHAPTER 3 PROBLEM STATEMENT A new computational method, namely, the material point method (MPM) was developed by Sulsky and others (1995,1996) at the University of New Mexico with support from the Sandia National Laboratories. Compared to FEM, MPM has the following advantages: (1) MPM is able to handle large deformation in a more natural manner so that mesh lockup present in FEM is avoided; (2) MPM can easily couple with molecular dynamics (MD) simulations because of the use of material points (similar to atoms used in MD) instead of arbitrary sized elements in FEM; (3) Parallel computation is more straightforward in MPM because of the use of a grid structure that is consistent with parallel computing grids; and (4) Use of the background grid in MPM enables structured adaptive refinement for local interested region. As a result, MPM is a logical choice for coupling MD to continuum for multiscale simulations. However, the conventional MPM discretizes the material using a regular grid with uniform square/rectangular (2D) or cubic/brick (3D) cells. Such a mesh is inefficient in dealing with stress concentration issues and 58 cannot be refined to scale down the size to different levels (Tan and Nairn 2002). Consequently, a hierarchical mesh refinement technique in MPM simulations needs to be developed. This study concentrates on the development of multiscale modeling technique. A seamless MD/MPM coupling approach will be developed in which MD modeling will be used in the local area of interest where the variation is at an atomistic scale and MPM simulation in the rest of the material. In this investigation, multiscale modeling and simulations of material properties will be addressed considering tensile testing. This application is chosen because of the extensive use of this technique for determining a wide range of material properties, including modulus, yield strength, ultimate strength, and strain at fracture. MPM/MD simulations of uniaxial tension will be carried on a silicon specimen at different strain rates to illustrate this methodology. A hierarchical mesh refinement technique will be introduced in MPM to scale down from the continuum level into the atomistic level near the MD/MPM transition region where the MPM points overlap the MD atoms. A seamless coupling is necessary for complete compatibility of both strain and stress for MD and MPM regions. Further, elastoplastic constitutive material model will be implemented in MPM for continuum scale simulation of plasticity material behavior. The framework has been established by Sulsky et al. (1996). Elastic plastic MPM at continuum scale coupled with MD at atomistic scale will aid in simulation of nanocrystalline material behavior. 59 CHAPTER 4 COMPUTATIONAL METHODOLOGY 4.1. Introduction. In coupled simulations of material behavior at different length scales, molecular dynamics (MD) is the preferred technique at atomistic scale while material point method (MPM) was chosen for simulation at continuum scale. A seamless coupling between the two techniques was recognized because of the natural similarity between the two methods. For a given body under analysis, molecular dynamics was used at the localized regions of interest, while the rest of the workmaterial is modeled using continuum mechanics approach. The two length scales are coupled by introducing a handshake region in between them, termed as ‘transition region’. A methodology to scale down the locally interested continuum region, to handshake with the atoms in the transition region, was achieved using hierarchical mesh refinement technique. The theoretical background along with computational methodologies adopted for coupled simulations are described in the sections that follow. 60 4.2 Molecular Dynamics (MD) 4.2.1 Theoretical Background In molecular dynamics (MD), the Newtonian second law for the time evolution of a system is expressed by (refer Komanduri et al. 2001,2003) 2 2 mid ri / dt fi , (i = 1,…, N) (4.1) where N is the total number of atoms, mi is the mass of the ith atom with position vector ri, i i 1 2 N f V(r ,r ,...,r ) is the force acting on the ith atom due to interaction with other atoms in the system, and 1 2 N V(r ,r ,...,r ) is the interaction potential. The initial positions and velocities of the atoms, together with the interaction potential, define the whole set of thermodynamic, elastic and mechanical properties of the material. Set of 3N secondorder differential equations [Eq. (4.1)] is often solved by recasting it as a set of 6N firstorder Hamiltonian equations of motion, thus i i i 1 2 N dp / dt f V(r ,r ,...,r ) , i i i dr / dt p /m . (4.2) where pi is momenta. Under the given initial positions and momenta of the system, integration of Eq. (4.2) yields the total trajectory of the system. With the knowledge of the trajectories of all the atoms, one can calculate spatial and temporal distributions of energy, temperature, and pressure, as well as monitor the structural and phase changes in the system. 61 The atomistic potential for the interaction between SiSi atoms herein is a general form of Tersofftype threebody potential (Tersoff 1988,1989). The total energy as a function of atomic coordinates is given by i j ij W 2 1 E , (4.3) where i and j label the atoms of the system and ij W is the bond energy of all the atomic bonds in the control volume represented by W f (r )[f (r ) b f (r )] ij c ij R ij ij A ij , (4.4) where ij r is the interatomic distance. The R f and A f functions are the repulsive and attractive forces between atoms i and j, respectively. These functions are expressed by f (r ) A exp( r ) R ij ij ij ij , (4.5) and f (r ) B exp( r ) A ij ij ij ij . (4.6) The f (r ) c ij function acts as a cutoff function that smoothly attenuates the interaction as ij r approaches the cutoff radius ij S . ij ij ij ij ij ij ij ij ij ij ij c ij 0 for r S ) for R r S S R r R cos( 2 1 2 1 1 for r R f (r ) (4.7) The ij b function is given by i i 2ni 1 n ij n ij ij i b (1 ) , (4.8) where 62 k i, j ij c ik ijk f (r )g( ) , (4.9) with 2 i ijk 2 i 2 i 2i 2i ijk d (h cos ) c d c g( ) 1 , (4.10) where ijk is the bond angle between bonds ij and ik. The parameters used in the Tersoff’s potentials for silicon, namely, A, B, R, S, , , , n, c, d and h are given in table 1 (Tersoff 1988,1989). Table 1. Material parameters in Tersoff’s potential for silicon (Tersoff 1988,1989) Parameter A (eV) B (eV) (Å1) (Å1) n c d h R (Å) S (Å) Silicon 1.83 103 4.71 102 2.48 1.73 1.1 105 0.79 1.0 105 16.22 0.6 2.7 3.0 4.2.2. Strain and Stress In MD simulations, stress and strain relation is usually determined from tensile models (Frankland et al. 2003, Chang and Fang 2003). In each increment of applied deformation, a nominal strain can be determined by the nominal length change normalized by initial length for the entire MD model in the loading direction. In thermodynamics, the stress in a solid is defined as the partial derivative of the internal energy with respect to the strain per unit volume. In continuum mechanics, the components of the stress tensor for a hyperelastic material are given by 63 ij S ij E V 1 , (4.11) where V is the volume of the solid, E is the total internal energy, ij is the components of the strain tensor, and S is entropy. In special cases, if the internal energy is the strain energy, the Hooke’s law is the result of Eq. (4.11). In molecular dynamics, the total energy is the summation of the energy of individual atoms, E , and is expressed as 1 2 E T U M ( ) ( ) 2 v r , (4.12) where T is the kinetic energy,U the potential energy, M the mass, v the magnitude of its velocity, and (r) the potential energy at the atom location r. Thus the stress components for a given atom are written as (M v v F r ) V 1 ij i j i j , (4.13) where V is the atomic volume of atom , i v and j v are components of velocity. i F are the components of force between atoms and from the derivative of the potential energy, and j r are the components of the distance between atoms and . The nominal stress is computed as the average of the atomic stresses for the volume of the model. Thus, the nominal stress components of each model are expressed by (M v v F r ) V 1 ij i j i j , (4.14) where V is the volume of the MD model and equals to V . 64 4.3. The Material Point Method (MPM) 4.3.1 Theoretical Background In the MPM (Sulsky et al. 1995, 1996) a material continuum is discretized into a finite collection of N material points, as shown in Fig. 14. Each material point is assigned a mass (mp, p=1, …, N) consistent with the material density and volume of the point, and all other variables, such as position, acceleration, velocity, strain, and stress. The motion equations are solved at the points for each time step of the analysis using a background computational grid. Information is transferred from the material points to grid nodes by using standard finite element shape functions, n i i i 1 ( ) N ( ) ! x ! x , (4.15) where n is the number of nodes in the grid and subscript i refers to the nodal values of !(x) containing x . Figure 14. MPM grid cells with material points. 65 The Newtonian equations of motion for material points are given by "d2u / dt2 " s "fb , (4.16) where the mass density N p p p 1 m ( ) " # x x , p m is material point mass, p x is material point location, #(xxp) is unit step function, u is the displacement, the acceleration a d2u / dt2 , b f is the specific body force, and s is a specific stress tensor, which is the stress tensor divided by the mass density. The material constitutive relation can be represented as s C(s) and 1 T ( ) ( ) 2 $& u u %' , (4.17) where C(s) is the fourthorder specific elasticity tensor, is the strain tensor. This paper will consider only linear elastic materials at the continuum level. The specific stress components are defined as s / " ij ij (4.18) Thus, Eq. (4.16) can be rewritten as i i s ij, j b a (4.19) Constitutive equations for a linear elastic material are given by k,l s ijkl s ij E v , (4.20a) with E E / " ijkl s ijkl (4.20b) where i v are the velocity components and ijkl E are components of the elasticity tensor. 66 Based on Eqs. (4.16)(4.20), the governing equations in weak form are ( ) () () () " ) " ) " ) T wT dS wb d wa d w d , j s ij i i i (4.21a) and ( E v ) wd 0 k,l s ijkl s ij ( " ) ) (4.21b) where w is the test function which is an arbitrary spatial function. Boundary conditions are u g i on u ) (on displacement boundary) and j s ij i T " n on T ) (on traction boundary), where j n are components of the unit outward normal vector at the boundary surface and u T ) ) ) . It is assumed that w 0 on u ) . A material continuum is divided into a finite collection of discrete infinitesimal regions p ) ( p p 1,...,N ) called material points. Each material point is assigned a mass p m in p ) , where () " ) p m (x)d p and p ) ) . Mass density can then be approximated as a sum of material point masses using a dirac function " # Np p 1 p p (x) m (x x ) (4.22) Applying Eq. (4.22) in Eq. (4.21a, b) converts integrals into summations of quantities evaluated at material points ( ) T p p Np p 1 (p) , j ) p ( , sp ij N p 1 (p) i (p) p N p 1 ) p (i (p) i p wT dS m w b m w a m w (4.23a) and 67 (p) k,l s ijkl ) p ( , sij E v . (4.23b) All variables represented by !(x) , such as the coordinates, displacements, velocities, and accelerations need to be transferred between grid nodes and material points using the shape functions N(x), ! ! N n 1 (x) (n)N(n) (x) . (4.24) where N is the number of nodes in the background grid and !(n) refers to the nth nodal values of !(x) . The numerical solution will be obtained at discrete set of times, t k , where k 1,...,K. Applying the interpolation scheme in Eq. (4.24) to the weak form of momentum equations in Eq. (4.23a) gives * * Np p 1 (p)(n) , j ) p ( k , sp ij n n' 1 N n 1 ) n ( k ) ' n ( ki k(nn') N n 1 wk(n) m a w m N * * N n 1 N n 1 N p 1 k(p) (p)(n) p i k(n) k(n) i k(n) p N b m w Tˆ w (4.25) Applying the interpolation scheme for acceleration into the weak form of Newton’s equation gives n int ext ij j i i j 1 m a f f , (4.26) where the mass matrix is given by N ij p i p j p p 1 m m N ( )N ( ) x x , (4.27) the internal and the external forces on the node i are 68 N int s i p p i p p 1 f m (x ) N (x ) and N ext i i p b p i p p 1 f ˆ m f (x )N (x ) + , (4.28) where i ˆ+ is the surface traction associated with grid point i. The grid point accelerations obtained from Eq. (4.26) are then used to update the position, velocity, stress, strain, and temperature of the material points. 4.3.2. Hierarchical Mesh Refinement (a) (b) Figure 15. Hierarchical mesh refinement in 2D MPM. The conventional MPM method uses a regular grid cell mesh with equally spaced nodes. Adaptive mesh refinement in areas with stress concentration is not possible using such a mesh. To accommodate multiple level meshes, a hierarchical mesh refinement technique is presented herein. For twodimensional 1 2 4 3 5 ,  1 2 4 3 ,  1 2 3 4 5 6 7 8 9 10 Level I Level II Level III 69 situation, the hierarchical mesh refinement with three levels using fournode regular cells (1, 4, 710) and fivenode transition cells (2, 3, 5, 6) are shown in Fig. 15 (a). For the regular cell with four nodes in Fig. 15 (b), the shape functions used in Eq. (4.15) are expressed as (Zienkiewicz et al. 1989,1991), (1 )(1 ) 4 1 N1 ,  , (1 )(1 ) 4 1 N2 ,  , (1 )(1 ) 4 1 N3 ,  , (1 )(1 ) 4 1 N4 ,  . (4.21) where (,, ) are natural coordinates, as shown in Fig. 15 (b). Shape functions for transition cells with five nodes are given by 1 5 N 2 1 (1 )(1 ) 4 1 N ,  , 2 5 N 2 1 (1 )(1 ) 4 1 N ,  , (1 )(1 ) 4 1 N3 ,  , (1 )(1 ) 4 1 N4 ,  , (1 )(1 ) 2 1 N 2 5 ,  . (4.22) Similarly, for threedimensional situation, a hierarchical refinement mesh with two levels is shown in Fig. 16 (a) using regular cells with eight nodes and transition cells with nine (cells 2, 4 and 6) and thirteen nodes (cells 1, 3 and 5). The shape functions for regular cells with eight nodes in Fig. 16 (b) are given by (Zienkiewicz et al. 1989, 1991) (1 )(1 )(1 ) 8 1 N i 0 0 0 ,  , 0 i , ,, , 0 i   , 0 i (i = 1,…,8), (4.23) i( , , ) i i i ,  : 1 (1, 1, 1); 2 (1, 1, 1); 3 (1, 1, 1); 4(1, 1, 1); 5(1, 1, 1); 6(1, 1, 1); 7(1, 1, 1); 8(1, 1, 1), and (,, , ) are natural coordinates as shown in Fig. 16. The shape functions for transition cells from nine to thirteen nodes in Fig. 16 (b) are expressed by 70 (1 )(1 )(1 ) 8 1 Nˆ i i i i ,,  (i = 1,…,8), 1 1 Nˆ N , 2 2 10 11 13 N 4 1 N 2 1 N 2 1 Nˆ N , 3 3 11 12 13 N 4 1 N 2 1 N 2 1 Nˆ N , 4 4 Nˆ N , 5 5 Nˆ N , 6 6 9 10 13 N 4 1 N 2 1 N 2 1 Nˆ N , 7 7 9 12 13 N 4 1 N 2 1 N 2 1 Nˆ N , 8 8 Nˆ N , (1 , )(1 )(1 ). / 4 1 N 2 9 , (1 ,)(1 )(1 ). ,/ 4 1 N 2 10 , (1 , )(1 )(1 ). / 4 1 N 2 11 , (1 ,)(1 )(1 ).,/ 4 1 N 2 12 , (1 , )(1 ).1 / 2 1 N 2 2 13 . (4.24) In 3D MPM computational formulation using the hierarchical mesh with multiple levels, Eq. (4.23) is used for regular cells at each individual level while Eq. (4.24) is employed for the transition cells between neighboring two levels. 71 (a) Regular cell Transition cell (b) Figure 16. 3D MPM regular and transition cells. 1 2 3 5 6 7 11 9 10 12 13  , 4 8 1 2 3 5 6 7  , 4 8 Level I I 6 5 4 Level I 1 2 3 72 4.4. Seamless Coupling between MD/MPM Simulation of dynamic problems using both MD and MPM requires integration of the equations of motion [Eqs. (4.1) and (4.16)]. The format of these two equations is similar, leading to a natural coupling of the two types of simulations. There are essentially three regions in coupled MD/MPM simulations. They are MD, MPM, and their transition region. The MD region is modeled using the approach described in Section 4.2. In MD region, atoms are initially placed at diamond lattice points. The MPM region is modeled using the method described in Section 4.3. In the MPM region, eight material points are equally positioned in each cubic cell. The coupling of the two descriptions of the media is brought about using a transition region where both MPM points and MD atoms are initially overlapped and coincide with the atomic lattice. The width of the transition region is equal to the cutoff distance of the interaction potential used in the MD region. This provides a complete set of neighbors within the interaction range for all atoms in the MD region. Atoms/points that belong to the transition region not only interact via the interaction potential with the MD region but also are part of points in the MPM region. The positions and velocities of these atoms/points in the transition region (both on the MD side and on the MPM side) must be consistent with each other. The atom velocities in the transition region due to the interaction with the MD region furnish the loading boundary condition on overlapped points for the MPM simulation. To avoid a density mismatch at the boundary, the mass at each point for MPM in the transition region is set equal to the mass of each MD atom. 73 Using the above method to deal with the transition region, the displacement as well as strain fields are compatible, while stress fields possibly are not compatible. To ensure complete compatibility of both strain and stress for seamless coupling between MD and MPM, we use the nonlocal elasticity theory (Abraham et al. 1998) in this study. The elastic energy is expanded in terms of the strain in the Taylor series form as [ E /( )] ... 6 1 [ E/( )] 2 1 E( ) E(0) ( E / ) ij kl mn 0 ij kl mn 3 ij kl 0 ij kl 2 ij 0 ij (4.25) Table 2. Material properties of single crystal silicon (Runyan 1999) Material properties Lattice constant (Å) Density (g cm3) Young’s modulus (GPa) C11 (GPa) C12 (GPa) C44 (GPa) Poisson’s ratio Bulk modulus (GPa) Fracture stress (MPa) Silicon 5.43 2.33 130 165.78 63.94 79.57 0.28 100 62 Therefore, the material elastic constants are required to satisfy the following conditions: ij ij 0 c ( E / ) , (firstorder elastic constants) ij kl 0 2 ijkl c [ E /( )] , (secondorder elastic constants) ij kl mn 0 3 ijklmn c [ E /( )] . (thirdorder elastic constants) 74 Using the natural state as the reference configuration, the firstorder elastic constants are zero and only the second and higherorder elastic constants need to be assigned to the continuum. As linear elasticity theory is used in this investigation, the secondorder elastic constants of the lattice and continuum regions require matching. Material properties of single crystal silicon are listed in table 2 (Runyan 1999). 75 CHAPTER 5 NUMERICAL IMPLEMENTATION OF SEAMLESS COUPLING BETWEEN MPM AND MD In this work, a tension model is used to conduct coupled MD/MPM simulations. Using the seamless coupling algorithms described in Chapter 4, a coupled MD/MPM code has been developed for multiscale simulations. The computational scheme or the flow chart used for the coupled MD and MPM simulation of silicon under tension is shown in Fig. 17. To generate the tension model, initial position data need to be assigned first to all atoms and points in both MD and MPM regions, which include atoms/points in the MD/MPM transition regions. Loading and boundary conditions are also prescribed for the entire system. Based on the given temperature, MD computation is first initiated so that the initial velocities for all atoms can be determined. The initial velocities will provide boundary conditions for the MPM region through the transition region. Then, MPM computation can be started and MPM results in the transition region provide boundary conditions for the MD region. Thus, MD computation can be iterated and MD results in the transition region provide new boundary conditions in MPM computation for the next time step. This iteration process continues until the completion of the simulation. For the purpose of illustrating the method, a 76 schematic of simple tension model, which consists of one MD region, two transition regions, and two MPM regions with three hierarchical mesh refinement levels, is shown in Fig. 18. Figure 17. Coupled MD and MPM computational flowchart. ti+1 = ti + 0t ( i = 0, 1, 2,….., n1) if ti+1 = tn Input data 1. Initial position r0 for atoms in MD region, atoms/points in transition regions, and points in MPM regions 2. Loading and boundary condition for entire system v 3. Temperature T for MD. Initial velocities vi0 for atoms in MD region and transition regions t = ti ( i = 0, 1, 2,….., n), 0t MPM computation ri+1, vi+1 for points/atoms in transition regions ri+1, vi+1 for points in MPM regions MD computation ri+1, vi+1 for atoms/points in transition regions ri+1, vi+1 for atoms in MD region END 77 Figure 18. Schematic of tensile model for coupled MD/MPM simulations. The next item of interest is coupling details in the transition region. For the MPM cells in the MD/MPM transition regions (shown in Fig. 18), MPM points are arranged following MD atoms at the atomic lattice sites, which have the interaction with the MD region. In other words, MD atoms and MPM points in the transition region are in onetoone correspondence and are initially overlapped at the same locations, as shown in Fig. 19 (a). It may be noted that MD time step increment is normally much smaller than that in MPM, so that MD and MPM computations have to be matched well in coupling simulations. We use adaptive time step increments in the transition region wherein the MD time step increment and MPM time step increment are the same. For continuum region, the time step increment gradually increases for the area farther away from the atomic region. In this simple model, we chose the same time step increment for both MD and MPM computations. For time integration, the time step increment is 0t = 1.018 fs MD/MPM transition region Atoms Points MPM MD region MPM 78 (1015 s). Details of the steps involved in coupled MD/MPM simulations are given below: (a) (b) (c) Figure 19 (a). MD atoms and MPM points initially overlapped in the transition region. (b). MPM points moving under velocity. (c). MD atoms moving same positions as MPM points. Black: MD atoms Grey: MPM points Dash line: MPM grid Step 1: Initial positions 0 r for all atoms in the MD region and all points in the MPM regions, which include atoms/points in the MD/MPM transition regions, are given according to the dimensions of the tension model and the crystal structure of the material. At the absolute temperature T, MD computation is initiated and initial velocities i0 v for all MD atoms are calculated. The initial velocities at atoms are applied to the overlapped MPM points in the transition regions as a boundary condition for MPM regions, as shown in Fig. 19 (b). Step 2: MPM computation is conducted at the time step t0 and new positions 1 r and velocities 1 v are obtained for all MPM points. The new positions r0 vi0 r1 r1 r0 v1’ 79 at points are applied to the overlapped MD atoms in the transition regions as a boundary condition for MD simulation as shown in Fig. 19 (c). Step 3: MD computation is performed at the time step 0 t and the updated positions r1 ' and velocities 1 v ' are obtained for all atoms. The updated velocities at atoms are applied to the overlapped MPM points in the transition regions as a new boundary condition for MPM simulations at the next time step 1 t . Step 4: Repeat the same algorithm to steps 2 and 3 in turns until the entire computation for the tension model is implemented by means of coupled MD and MPM approach. During coupling of MD with MPM simulations, MPM always provides the updated positions for MD atoms in the transition regions, while MD in turn provides the updated velocities for MPM points in the transition regions. It should be noted that following this algorithm from step 1 to step 4, both atoms for MD and points for MPM in the transition region naturally maintain displacement compatibility so that strains are compatible as well. However, the stress compatibility cannot be ensured automatically. Additionally, mass mismatch exists for both MD atoms and MPM points in the transition region, since mass determined by the mass density and the occupied volume is assigned to each point in the MPM region. To ensure complete compatibility of both strain and stress for seamless coupling between MD and MPM, nonlocal elasticity theory is 80 used in this study. Based on this theory, for homogenous, linear elastic materials, only the secondorder elastic constants require matching between atomistic and continuum regions. In this investigation, we take elastic constants for continuum region from the lattice theory and use the same mass for both MD atoms and MPM points in the transition regions. 81 CHAPTER 6 MULTISCALE SIMULATIONS OF TENSILE TESTING 6.1. Introduction In this investigation, multiscale modeling and simulations of material properties are addressed considering tensile testing. This application is chosen because of the extensive use of this technique for determining a wide range of material properties, including modulus, yield strength, ultimate strength, and strain at fracture. In this work, a tension model is used to conduct coupled MD/MPM simulations. Using the seamless coupling approach, multiscale simulations of tensile testing are carried on single crystal silicon workmaterial. Results are presented in terms of stress  strain relationships at several strain rates, as well as the rate dependence of uniaxial material properties. The simulation details and the results obtained are given in the sections that follow. 82 6.2. Coupled MD/MPM tensile testing model. Fig. 20 is 3D animation of coupled MD/MPM model at the initial stage of tensile testing showing the MD region, MD/MPM transition regions, and MPM regions. There are three levels in the MPM region. The dimensions of this model are 32.6 Å 32.6 Å 222.8 Å. There are 2,901 atoms in the MD region and 9,280 material points in the MPM regions, which include atoms/points in the MD/MPM transition regions. In the transition regions, MD atoms and MPM points are initially overlapped in onetoone correspondence and coincide with the atomic lattice. In the MD region, the crystal is set up with a cube orientation [001] and the tensile loading is applied along the [001] direction under controlled velocity. The initial temperature is 293 K for the MD simulations. The elastic constants, the Young’s modulus, the Poisson’s ratio, and other relevant parameters, given in Table 2, are used in the MPM regions. The coupling simulations on the tension model are conducted for the silicon workmaterial at different strain rates, namely, 44.88 ns1, 53.86 ns1, 62.84 ns1, 67.33 ns1, 80.79 ns1, and 89.76 ns1. It may be noted that very high strain rates are used of necessity to minimize the computational time for coupling simulations. For the coupling simulation at a strain rate of 44.88 ns1, the total simulation time was 4 ps (1012 s) and the time step increment used was 1.018 fs. It took ~ 3 hrs of computational time using a PC with a 2.8 GHz processor. 83 Figure 20. Animation of coupled MD/MPM model showing initial stage of tensile testing. MD Atoms MPM Points MPM Points MD region Level I MD/MPM transition region Level II Level III MPM region Level III Level II Level I MPM region 84 6.3. Deformation and failure mechanism under tensile loading Fig. 21 shows snapshots of coupling simulation for the entire model at various stages of tensile testing at a strain rate of 44.88 ns1 and Fig. 22 shows closeup views of Fig. 21. The animations of coupled MD/MPM simulation show that the tensile specimen experiences with increasing elongation in elasticity, generation of dislocations and plasticity by slip, void formation and propagation, formation of amorphous structure, necking, and final fracture in the MD region while only elastic deformation is present in the MPM regions at either end. In the initial stage, there is uniform stretching with elastic deformation and no dislocations occurring in the crystal lattice, as shown in Fig. 23 (b). The limit of elasticity is reached when the strain reaches approximately 0.1, as shown in Fig. 27. With further increase in deformation, the dislocations are seen to be generated locally with their densities gradually increasing in the plastic deformation region. Slip takes place along [111] at about 45o to the loading direction as shown in Fig. 23 (c). Plasticity by slip lasts until the formation of voids and semiamorphous structure, as shown in Fig. 23 (d). Ultimate tensile stress is reached at a strain of ~ 0.3, as shown in Fig. 27. Further increase in elongation results in local interactions or rearrangements of silicon atoms, resulting in somewhat amorphous structure [see Figs. 23 (e) and (f)]. The deformation is then concentrated in a narrow region resulting in necking, as shown in Fig. 21 (d). Finally, separation of the tensile specimen occurs as shown in Fig. 21 (e). 85 (a) (b) (c) (d) (e) Figure 21. Snapshots of MD/MPM coupling simulation at strain rate 44.88 ns1. 86 ( a ) ( b ) ( c ) ( d ) ( e ) Figure 22. Snapshots of MD/MPM coupling simulation at strain rate 44.88 ns1 (magnified view). 87 Figs. 22 (a) to (e) are 3D animations of MD simulation at higher magnification (of Fig. 21) showing different stages of tensile testing at a strain rate of 44.88 ns1. Slip taking place at 45o to the loading direction can be clearly seen [Figs. 22 (b) to (c)]. The formation of amorphous structure in the central region and void formation can also be clearly seen in Fig. 22 (d). Fig. 22 (e) shows concentration of amorphous structure in the central region followed by necking and separation. This corresponds to the steepest slope in Fig. 27 leading to final rupture of the sample. Figs. 23 to 27 show animations in the midsection at different stages of tensile testing at different strain rates from 44.88 ns1 to 89.76 ns1. Initially the deformation is predominantly elastic. This corresponds to steep elastic slope in Fig. 27. This is followed by slip on the preferred [111] direction. Further increase in strain leads to the formation of inhomogeneous structure leading to amorphization in the entire necked region [see atoms marked by circles in Fig. 23 (d)]. A similar amorphization can be seen in Figs. 24 (d), 25 (d) and 26 (d). Ultimately, separation of the sample takes place with drastic reduction in stress. With increase in strain rate, the crystallinity of the silicon is maintained up to larger strains. Also, with increase in strain rate, the strain at fracture increases. It may be noted that necking that occurs before final failure for silicon is significantly less as compared to ductile materials, such as copper, aluminum (Komanduri et al. 2001). The catastrophic failure leads to a sudden drop in the stressstrain curves as shown in Fig. 27. 88 (a) (b) (c) (d) (e) (f) Figure 23. Snapshots of MD/MPM coupling simulation at strain rate 44.88 ns1. 89 (a) (b) (c) (d) (e) (f) Figure 24. Snapshots of MD/MPM coupling simulation at strain rate 67.33 ns1. 90 (a) (b) (c) (d) (e) (f) Figure 25. Snapshots of MD/MPM coupling simulation at strain rate 80.79 ns1. 91 (a) (b) (c) (d) (e) (f) Figure 26. Snapshots of MD/MPM coupling simulation at strain rate 89.76 ns1. 92 6.4 Effects of strain rate on material properties The effects of strain rate on material properties have also been investigated for the coupled MD/MPM simulations. Under different strain rates, slip essentially takes place at about 450 with respect to the loading direction followed by the formation of amorphous structure and generation of voids in the necked region and their coalescence, leading to final failure. The engineering stressstrain curves at different strain rates are shown in Fig. 27. They show that at different strain rates material behaves differently at different stages including elastic and plastic deformation, and final rupture. Thus, the strain rate affects material properties significantly. However, before the strain reaches ~ 0.1, there is no dislocation generation and the initial part of engineering stressstrain curve is essentially elastic at different strain rates in the range from 44.88 ns1 to 89.76 ns1. From the slope of the stressstrain curves, the Young’s modulus is determined as ~ 121 GPa, which is very close to that obtained from the lattice theory (Runyan, 1999). As expected, the elastic modulus is nearly constant and is not affected by the strain rate. The curve in Fig. 28 (a) shows the effect of the strain rate on the tensile strength for silicon. The tensile strength increases from ~ 15 to ~ 18 GPa with increase in the strain rate from 44.88 ns1 to 89.76 ns1. Fig. 28 (b) shows the variation of the ultimate strain with strain rate. It shows that the ultimate strain increases rapidly up to ~ 0.54 as the strain rate increases from 44.88 ns1 to 67.33 ns1, and then increases gradually up to ~ 0.57 as the strain rate increases from 67.33 ns1 to 89.76 ns1. The ultimate strain vs. the strain rate curve 93 (GPa) 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0 5 10 15 20 53.86 ns 67.33 ns1 1 44.88 ns1 62.84 ns1 80.79 ns1 89.76 ns1 Figure 27. Engineering stressstrain curves from coupling MD/MPM simulations for silicon. changes slope at a strain rate of ~ 67.33 ns1. Observation of the animations indicates that defect initiation and coalescence are somewhat different across the strain rate. At lower strain rates, i.e. below 67.33 ns1, [for example at 44.88 ns1, as shown in Fig. 23(d)], voids are formed in the central region of the tensile specimen. These voids coalesce to larger defects resulting in necking and ultimate failure. While at higher strain rates (67.33 ns1 or higher), [for example at 89.76 ns1, as shown in Figs. 26(d) and (e)], defects are formed near the surface at about ¼ and ¾ of the MD region (see regions marked by circles in Fig. 26). These defects grow in size and cover the entire area, leading to decreased load bearing capacities at these locations, resulting in some bulging, as shown in Figs. 26 (d)  (e). Compared with failure at the lower strain rates 94 40 50 60 70 80 90 12 14 16 18 20 (ns ) 1 u (GPa) (a) Tensile strength vs. strain rate. 40 50 60 70 80 90 0.3 0.4 0.5 0.6 0.7 (ns ) 1 u (b) Ultimate strain vs. strain rate. Figure 28. Effects of strain rate on material properties. 95 [e.g., Fig. 23(d)], there are more defect sites at higher strain rate [Fig. 26(e)]. Once the defects at multiple sites coalesce, necking occurs, leading to a quick rupture of the material. The change in the slope of the ultimate failure strain shown in Fig. 28 (b) is most likely the result of the change in defect formation mechanism from relatively low to high strain rates. However, it is not clear at this stage the underlying reason for a change in the defect formation mechanism across the strain rates. Further investigation is needed to explore this feature. 96 CHAPTER 7 ELASTOPLASTIC ANALYSIS USING MATERIAL POINT METHOD 7.1.1 Metal plasticity models In processing of metallic materials, elasticplastic constitutive material behavior is involved. Numerical analysis for this class of materials is an important aspect. A wide range of processes, such as metal cutting, metal forming involve large deformation and plastic flow of the material. The material point method (MPM) has been recently found useful for such applications (Sulsky et al. 1994). In plasticity, numerous constitutive models exist and the main choices are between rate independent and rate dependent plasticity, between the von Mises and Hill’s yield surface. Further, we have choice within rate independent model as either isotropic hardening or kinematic hardening. Rate independent model is adequate for modeling the response of metals at low strain rates, and at operating temperatures much lower than the melting temperature. For metals, yielding is independent of the hydrostatic stress. Experiments advocate that von Mises yield criterion provides somewhat better fit to experimental data than Tresca yield criterion. Hence von Mises yield surface is used in this formulation to define isotropic yielding. A state of stress on the von 97 Mises yield surface causes plastic deformation and yielding. In case strain hardening occurs, the radius of the yield surface increases. In the case where polycrystalline metals obey isotropic plastic hardening law, the radius of the yield surface changes in all directions. In dynamic problems such as tension and nanoindentation, which involve finite strains, as well as in manufacturing processes, such as metal forming and metal cutting, large plastic strains are present and plastic strain does not continuously reverse direction sharply. In such applications, isotropic hardening theory is usually considered. Associated plastic flow rule is used, which means that the plastic deformation increment has the same direction as the outward normal to the yield surface. It should be noted that, there are also some limitations posed by the model, namely, the model does not predict behaviour under cyclic loading correctly. Since the incremental plasticity model is based on the assumption of small displacements and small strains, it will not predict plastic strains correctly if the principal axes of stress rotate significantly during deformation. The general equations for the plasticity model using von Mises yield surface, isotropic hardening and associative flow rule, are summarized in the following section. Two examples, one, a three  dimensional tensile bar under impact load, and the other Taylor Impact bar problem are used to validate the code. For the later, results are presented to demonstrate axisymmetric formulation of the material point method (Sulsky and Schreyer, 1996). In computational mechanics finite element method (FEM) is being used for a wide range of problems with satisfactory results. Realizing the similarity in 98 basic dynamic equations that are solved by both MPM and FEM and that too using explicit formulation, FEM results are chosen for comparison and validation of MPM elastoplastic constitutive material model. 7.1.2 Summary of rate independent multiaxial plasticity model with von Mises yield surface and isotropic hardening. 1 Using Incremental approach and additive strain decomposition: pi j e ij ij , dt d dt d E dt d plm mn n ijmn ij , where, ij : Total incremental 2nd order strain tensor. t : time. e ij : Elastic part of the ij . pi j : Plastic part of the ij . Eijmn: 4th order elasticity tensor. ij : 2nd order stress tensor. 1 von Mises yield criterion, f( ij ) : 0 p ij SijSij H 2 3 f( ) , ij ij kk ij 3 1 S , where, ij S : Deviotoric part of ij . kk : 1st invariant of ij . 99 p : Equivalent plastic strain. p H : Radius of the yield surface. ij : 2nd order unit tensor. Also, pi j pi j p , where, pi j : Components of plastic strain tensor. 1 Isotropic strain hardening law: Linear hardening: . p / H = 0 + c p . Power law hardening: . p / H = 0 + c . /p m . where, 0 : Virgin yield stress in uniaxial stressstrain behavior. m: material constant. c: Slope of yield stress versus plastic strain curve, positive and negative values of which corresponds to hardening or softening, respectively, of the material. 1 Associative plastic flow rule: mn mn ij ij pi j S S S dt d d df dt d dt d 2 2 , where, dt: time increment. 2 : Measure of the plastic strain path length, and is a nonnegative parameter. 1 Loading / unloading condition in Kuhn  Tucker complimentary form: Find trial stress increment, 100 0 trial = E 0 and initialize 0 0 . If f( ) 0 and 0 d df 0 3' % 4& $ , then instantaneous elastic process, else If f( ) = 0 and 0 5 0 3' % 4& $ , then instantaneous plastic process. The loading / unloading condition above can be solely determined using trial state only, which is ultimately implemented in return mapping algorithm. 7.1.3 Radial return mapping algorithm for threedimensional nonlinear isotropic hardening plasticity. 7.1.3.1 Main algorithm. 1) Compute trial elastic stress and deviatoric stress, assuming all the strain increment is elastic. 6tr. /7 I 3 1 n 1 n 1 n 1 e , . / pn n 1 trial n 1 s 2 e e , where, I : Unit tensor; n 1 e : deviatoric part of total strain. n 1 : strain at current time; trial1 n s : trial deviatoric stress; K: Material bulk modulus; : shear modulus; p n e : plastic part of total strain. Subscripts (n+1) and (n) indicate current and previous time step values respectively. 101 2) Check yield condition. . p / n trial1 n trial n 1 H 3 2 f s IF f trial 0 n 1 then all the stress increment is elastic so update total stress ( ) and exit. . / trial1 n 1 n 1 n Ktr I s , exit. ELSE plas 



A 

B 

C 

D 

E 

F 

I 

J 

K 

L 

O 

P 

R 

S 

T 

U 

V 

W 


