

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


MULTISCALE SIMULATION USING THE GENERALIZED INTERPOLATION MATERIAL POINT METHOD, DISCRETE DISLOCATIONS AND MOLECULAR DYNAMICS By JIN MA BACHELOR OF SCIENCE Wuhan University of Technology Wuhan, China June, 2000 MASTER OF SCIENCE Oklahoma State University Stillwater, Oklahoma August, 2002 Submitted to the Faculty of the Graduate College of Oklahoma State University in partial fulfillment of the requirements for the Degree of DOCTOR OF PHILOSOPHY May, 2006 ii MULTISCALE SIMULATION USING THE GENERALIZED INTERPOLATION MATERIAL POINT METHOD, DISCRETE DISLOCATIONS AND MOLECULAR DYNAMICS Dean of the Graduate College Thesis Advisor Dissertation Approved: Hongbing Lu, Ranga Komanduri Jay C. Hanan Lionel M. Raff Demirkan Coker A. Gordon Emslie iii © Copyright By JIN MA May 05, 2006 iv Acknowledgements I wish to express my most sincere appreciation to my thesis advisor, Dr. Hongbing Lu, for his intelligent supervision, constructive guidance, inspiration and friendship. Dr. Lu has also given me very generous financial support. My sincere appreciation extends to other committee members, Dr. Ranga Komanduri, Dr. Demir Coker and Dr. Lionel Raff, whose guidance, assistance, encouragement, and friendship are also invaluable. I would like to thank Dr. Bo Wang, Dr. Jay Hanan, Dr. Samit Roy, Dr. Martin Hagan, Dr. Afshin J. Ghajar, Dr. C. Eric Price and Dr. J. Keith Good for their long time assistance and guidance. Moreover, I wish to express my sincere gratitude to those provided suggestions and assistance for this study: Mr. Gang Huang, Ms. Haowen Yu, Mr. Nitin Daphalapurkar and Dr. Huiyang Luo. My research work was supported by a grant from the Air Force Office of Scientific Research (AFOSR) through a DEPSCoR grant (No. F496200310281). I would like to thank Dr. Craig S. Hartley, and Dr. Jamie Tiley, Program Managers for the Metallic Materials Program at AFOSR for the support of and interest in this work. Thanks are also due to Dr. Richard Hornung and Dr. Andy Wissink of the Lawrence Livermore National Laboratory for providing the code and assistance with the use of SAMRAI, Dr. Scott Bardenhagen for his introduction to the GIMP algorithm and valuable comments, and Dr. John A. Nairn for providing the 2D MPM code. I want to give my special appreciation to my wife, Yang Liu, for her precious suggestions v to my study, her encouragement at times of difficulty, her love and understanding throughout my study. Thanks also go to my parents and my friends for their support and encouragement. Finally, I would like to thank the School of Mechanical and Aerospace Engineering for supporting me during my M.S. and Ph.D. studies. vi Table of Contents Acknowledgements............................................................................................................ iv Chapter 1 Introduction ........................................................................................................ 1 Chapter 2 Multiscale Simulations Using Generalized Interpolation Material Point (GIMP) Method and SAMRAI Parallel Processing ......................................................................... 4 2.1 Introduction.........................................................................................................5 2.2 Generalized Interpolation Material Point (GIMP) Method ................................ 9 2.2.1 Contact algorithm in GIMP ...................................................................... 12 2.2.2 Numerical implementation........................................................................16 2.3 Parallel Computing Scheme Using GIMP with SAMRAI ............................... 17 2.3.1 SAMRAI...................................................................................................17 2.3.2 Spatial and temporal refinements.............................................................. 18 2.3.3 Domain decomposition .............................................................................22 2.4 Numerical Examples and Results ..................................................................... 25 2.5 Conclusions.......................................................................................................34 Chapter 3 Structured Mesh Refinement in GIMP Method for Simulation of Dynamic Problems ........................................................................................................................... 37 3.1 Introduction.......................................................................................................38 3.2 GIMP.................................................................................................................40 3.3 Structured Mesh Refinement ............................................................................44 3.4 Numerical Simulations......................................................................................47 3.4.1 Tracking particle deformation................................................................... 47 3.4.2 Indentation problem..................................................................................51 3.4.3 Verification of the displacement boundary condition............................... 54 3.4.4 Stress concentration problem.................................................................... 56 3.4.5 Static stress intensity factor ...................................................................... 58 3.5 Conclusions.......................................................................................................62 Chapter 4 Simulation of Deformation Mechanism and Failure of Bulk Metallic Glass (BMG) Foam Using the GIMP Method............................................................................ 64 4.1 Introduction.......................................................................................................65 4.2 Nanoindentation Testing of BMG Foam .......................................................... 70 4.3 Reconstruction of the BMG Foam and Simulation Convergence .................... 72 vii 4.3.1 Reconstruction of the BMG foam for GIMP ............................................ 72 4.3.2 Effect of grid cell size on simulation convergence................................... 74 4.4 Simulation of the BMG Foam in Compression ................................................ 76 4.4.1 Material failure..........................................................................................78 4.4.2 Stress flow in the foam ............................................................................. 79 4.4.3 The minimum size of a representative volume element (RVE)................ 81 4.4.4 Simulation of densification ....................................................................... 84 4.5 Conclusions.......................................................................................................89 Chapter 5 Multiscale Simulation Using GIMP Method and Molecular Dynamics .......... 91 5.1 Introduction.......................................................................................................92 5.2 GIMP and Refinement ...................................................................................... 96 5.3 MD Simulation and Atomistic Strain ............................................................... 99 5.3.1 Incremental atomic strain........................................................................100 5.3.2 Interpolation function..............................................................................101 5.3.3 Numerical validation...............................................................................104 5.4 Coupling of GIMP and MD Simulations ........................................................ 107 5.4.1 Coupling scheme.....................................................................................107 5.4.2 Parallel processing of the coupled model ............................................... 111 5.5 Numerical Results........................................................................................... 112 5.5.1 Multiscale simulation of mode I crack.................................................... 112 5.5.2 Multiscale simulation of mode II crack .................................................. 117 5.6 Conclusions..................................................................................................... 119 Chapter 6 Multiscale Simulation of Nanoindentation using GIMP, Dislocations Dynamics and Molecular Dynamics ................................................................................................ 122 6.1 Introduction..................................................................................................... 122 6.2 Coupling Scheme between GIMP, DD and MD............................................. 126 6.2.1 Coupling of GIMP and MD .................................................................... 126 6.2.2 Bridging the continuum and atomistic scales with DD .......................... 128 6.2.3 Parallel processing ..................................................................................130 6.2.4 Contact solver for MD ............................................................................ 133 6.3 Simulation of Nanoindentation ....................................................................... 133 6.4 Conclusions..................................................................................................... 140 Chapter 7 Summary and Future Work............................................................................ 141 7.1 Summary ......................................................................................................... 141 7.2 Future Work.................................................................................................... 143 7.2.1 Other multiscale problems ...................................................................... 143 7.2.2 Simulation of metal cutting using GIMP................................................ 144 viii Disclaimer on SAMRAI ................................................................................................. 149 References...................................................................................................................... 150 ix List of Figures Fig. 21: Illustration of early contact in multimesh ......................................................... 13 Fig. 22: Schematic of contact algorithm between the rigid indenter and the deformable workpiece......................................................................................................................... 14 Fig. 23: Material points in cells and the weighting function in 2D GIMP ...................... 16 Fig. 24: Illustration of a hierarchy of three nested grid levels of mesh refinement......... 18 Fig. 25: Two neighboring coarse and fine grid levels in 2D GIMP computations.......... 21 Fig. 26: Nested multilevel refinement and reduction in the number of cells with number of levels............................................................................................................................ 21 Fig. 27: A computational domain of two patches in one grid level................................. 23 Fig. 28: Flowchart showing advancement of grid levels recursively starting from the coarsest to the finest level in GIMP.................................................................................. 24 Fig. 29: Simulation results of tensile stress contours for a simple tension problem ....... 26 Fig. 210: Loading conditions for a simple indentation problem ..................................... 27 Fig. 211: GIMP and FEM results of normal stress variation in the Ydirection at different increments .......................................................................................................... 28 Fig. 212: Schematic of 2D indentation showing (a) three levels of refinement and (b) the indenter velocity history ................................................................................................... 29 Fig. 213: Comparison of normal stresses in Ydirection from FEM and GIMP............. 30 Fig. 214: Comparison of shear stresses from FEM and GIMP ....................................... 31 Fig. 215: Normal and shear stresses of GIMP simulations with a uniform cell size of 500 nm .............................................................................................................................. 32 Fig. 216: Indentation load versus depth curves from FEM and GIMP with different grid sizes.................................................................................................................................. 33 Fig. 217: Multiscale simulation of nanoindentation with five levels of refinement........ 34 Fig. 31: 2D representation of a particle and a grid cell ................................................... 44 Fig. 32: Refinement of structured grid with a refinement ratio of two ........................... 45 Fig. 33: Effect of influence zone on the weighting functions ......................................... 47 Fig. 34: Schematics showing overlaps and gaps that may occur when particle deformation is not tracked ................................................................................................ 48 Fig. 35: Simulation showing separation when the particle deformation is tracked by strain................................................................................................................................. 49 Fig. 36: Velocity field of a continuum region (the arrows represent both direction and magnitude) ........................................................................................................................ 50 x Fig. 37: GIMP results with tracking deformation of corners and their comparison with FEM ................................................................................................................................. 51 Fig. 38: Twodimensional indentation simulation........................................................... 52 Fig. 39: Comparison of the stress distributions at different levels of refinements in force indentation........................................................................................................................ 53 Fig. 310: Comparison of displacement history with FE.................................................. 54 Fig. 311: A simple tension problem with displacement boundary conditions ................ 55 Fig. 312: Comparisons of the displacement and stress.................................................... 56 Fig. 313: A plate with a circular hole subjected to tension ............................................. 56 Fig. 314: Normal stress in the Xdirection with p = 10 MPa .......................................... 57 Fig. 315: Normalized tangential stress along the circumference of the hole .................. 57 Fig. 316: Geometry of a double cantilever beam with a crack........................................ 58 Fig. 317: Material points and background grid around the crack tip .............................. 60 Fig. 318: Computed stress intensity factor as a function of time .................................... 61 Fig. 319: Comparison of the stress intensity factor with MPM results .......................... 61 Fig. 320: Stress distribution in the beam using three levels of refinements at t = 2.5 ms62 Fig. 41: Illustrations of the discretization scheme and intrinsic contact in GIMP .......... 68 Fig. 42: High resolution TEM images of (a) crystalline and (b) amorphous structures of a zirconium alloy (Hufnagel (2005)) ................................................................................ 69 Fig. 43: Average nanoindentation loaddisplacement curve with error bars for the BMG foam (Pd43Ni10Cu27P20)..................................................................................................... 71 Fig. 44: Schematic of Xray tomography system ............................................................ 72 Fig. 45: 3D microtomographic image of a BMG foam................................................... 73 Fig. 46: Stress distribution on a section simulated with the grid size of 3 and 4 times of the voxel size, respectively ............................................................................................... 75 Fig. 47: Stress distribution on a section simulated with the grid size to be 1 and 2 times of the voxel size. ............................................................................................................... 76 Fig. 48: Experimental setup and nominal stressstrain curve.......................................... 77 Fig. 49: Stress distribution in the overall GIMP model at t = 2.5 μs at location 1.......... 79 Fig. 410: Stress distribution on different sections at t = 2.5 μs. ...................................... 80 Fig. 411: Stress distribution on different sections at 8% strain. ...................................... 80 Fig. 412: Stressstrain curves of the BMG foam from experiments and simulations using second order interpolation................................................................................................. 82 Fig. 413: Stressstrain curves of the BMG foam from experiments, model prediction and simulations using linear interpolation............................................................................... 83 Fig. 414: Stressstrain curve from the densification simulation...................................... 85 Fig. 415: Comparison of initial and final shapes of the RVE.......................................... 86 Fig. 416: Maximum shear stress on section Y50 at different compressive strains ........ 87 Fig. 417: Distribution of z σ stress on section Y50 at different compressive strains .... 88 xi Fig. 418: Histograms of z ε at three compressive strains ................................................ 89 Fig. 51: Two neighboring coarse and fine grid levels in 2D GIMP computations.......... 98 Fig. 52: 2D representation of a particle and a grid cell ................................................. 103 Fig. 53: The generalized interpolation function in 2D .................................................. 104 Fig. 54: Two examples to calculate the atomistic strain: (a) Tension (b) Shear ........... 105 Fig. 55: Strain histories from different simulations....................................................... 106 Fig. 56: Comparison of the normal strain of a surface atom in tension ........................ 106 Fig. 57: Illustration of coupled GIMP and MD simulations. The circles represent atoms and squares (smaller than physical size) material points. The material points connect to each other without a gap to represent continuum. .......................................................... 108 Fig. 58: Flow chart of the coupling algorithm for one increment ................................. 110 Fig. 59: Illustration of three GIMP grid refinements and the domain decomposition for coupling.......................................................................................................................... 112 Fig. 510: Coupled GIMP/MD simulation of a 2D modeI edge crack. The dashed lines are the boundaries of the atomistic domain .................................................................... 112 Fig. 511: Stress distribution for P = 0.3 eV/Å3.............................................................. 114 Fig. 512: Simulation model with an edge crack with P = 0.3 eV/Å3 at t = 208 τ.......... 115 Fig. 513: Energies in the model with P = 0.15 eV/Å3 ................................................... 116 Fig. 514: Energy release rate P = 0.15 eV/Å3................................................................ 116 Fig. 515: Boundary conditions of the mode II crack problem....................................... 117 Fig. 516: Comparison of results from pure MD and coupled simulations at t = 11.2 τ 118 Fig. 517: Dislocation path at different simulation time................................................. 119 Fig. 61: Illustration of coupled GIMP and MD simulations. The circles represent atoms and squares (smaller than physical size) material points. The material points connect to each other without a gap to represent continuum ........................................................... 127 Fig. 62: Illustration of the domain decomposition and refinement for the coupling simulation of 2D indentation using GIMP, DD and MD................................................ 131 Fig. 63: Flow chart of the coupling algorithm for each increment................................ 132 Fig. 64: 2D modeling of an indentation problem with a cylindrical indenter ............... 134 Fig. 65: Dislocations and the loaddepth curve ............................................................. 134 Fig. 66: Comparison of loaddepth curves using different contact algorithms ............. 135 Fig. 67: Effect of the indentation velocity on the loaddepth curve .............................. 136 Fig. 68: Effect of the indenter size on the loaddepth curve.......................................... 137 Fig. 69: Loaddepth curves when the ratio between the workpiece size and the indenter radius is constant............................................................................................................. 138 Fig. 610: Comparison of the loaddepth curves between the coupling and purely MD simulations ...................................................................................................................... 139 Fig. 611: Slip bands, discrete dislocations and stress distributions in the model.......... 139 Fig. 71: Simulation snap shots of 2D orthogonal metal cutting using GIMP ............... 147 xii List of Tables Table 41: Nanoindentation results for the BMG foam. ................................................... 71 Table 42: Porosity of different volume elements at each location. ................................. 81 Table 71: JohnsonCook material parameters for copper (ABAQUS (2005)).............. 146 Table 72: Material properties for copper....................................................................... 146 1 Chapter 1 Introduction Multiscale simulation is a very important tool to reveal material behavior at different length and temporal scales through numerical simulations. Despite significant developments in materials simulation techniques, the goal of predicting reliably the properties and behaviors of new materials has not yet been achieved. This situation exists for several reasons that include a lack of full understanding of material behavior at different scales, absence of scaling laws, computational limitations, and difficulties associated with experimental measurements of material properties at micro and nano scales. For example, the information on the mechanical behavior of materials at nano level is not presently available as input to nanotechnology for the manufacturing of nanocomponents or microelectromechanical systems (MEMS). Scaling laws governing the mechanical behavior of materials from atomistic (nano), via mesoplastic (micro), to continuum (macro) scales are very important to numerous 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. Scaling laws are also important for 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 2 conventional indentation. Appropriate scaling laws may extend the extensive knowledge accumulated over time on material behavior at the macro (or continuum) level to the atomistic (or nano) level, via mesoplastic (micro) level. This work is focused primarily on developing coupling methodologies for continuum and atomistic simulations. At the continuum level, a new simulation method, the generalized interpolation material point (GIMP) method (Bardenhagen and Kober (2004)) is used. GIMP is developed based on the conventional material point method (MPM) (Sulsky et al. (1995)) with improved simulation stability. In order for GIMP to cover several length scales, and to couple MD simulations, refinement techniques must be developed. The thesis is divided into several chapters. Each chapter includes an introduction specific to that topic and a conclusion of findings. Chapter 2 describes the general theory of GIMP and multilevel refinement technique at the continuum level based on the Structured Adaptive Mesh Refinement Application Infrastructure (SAMRAI) (Hornung and Kohn (2002)). Parallel processing is implemented for GIMP based on SAMRAI and a contact algorithm is developed for the simulation of nanoindentation problems. Chapter 3 presents a structured mesh refinement algorithm for GIMP that is smoother than the method in Chapter 2. The new refinement technique is based on structured mesh and the refined mesh remains structured. The displacement boundary condition is explicitly introduced into the discretization of equation of motion in GIMP. Chapter 4 describes an application of the GIMP method for 3D simulations of bulk metallic glass (BMG) foams prepared by thermoplastic expansion. The material points are created based on the voxels in Xray tomography so that the microstructure of the foam is reconstructed exactly in GIMP. The simulation results of uniaxial compression 3 are compared with experimental ones. The densification of foam during compression is analyzed. Chapter 5 presents a coupling algorithm for GIMP and MD simulations. A new method to compute the atomic strain from MD simulation was introduced and validated. The coupling algorithm enables force, velocity, displacement, and energy compatibility at the handshaking region of the continuum and atomistic domains. Numerical simulations are performed to validate the coupling algorithm. The coupling algorithm Chapter 5 is extended to include discrete dislocations in Chapter 6. The discrete dislocation field is coupled with the continuum field using superposition of boundary conditions. The dislocations in the MD simulation can be detected and passed into the continuum region and vice versa. The summary and future works are presented in Chapter 7. The preliminary results of using GIMP to simulate metal cutting is presented in Chapter 7. This problem involves issues such as contact, large deformation, material failure, and adaptive mesh refinement. 4 Chapter 2 Multiscale Simulations Using Generalized Interpolation Material Point (GIMP) Method and SAMRAI Parallel Processing In the simulation of a wide range of mechanics problems including impact/contact/penetration and fracture, the material point method (MPM) (Sulsky, Zhou and Shreyer, (1995)) demonstrated its computational capabilities. To resolve alternating stress sign and instability problems associated with conventional MPM, Bardenhagen and Kober (2004) introduced recently the generalized interpolation material point (GIMP) method and implemented for onedimensional simulations. In this chapter GIMP has been extended to 2D and applied to simulate simple tension and indentation problems. For simulations spanning multiple length scales, based on the continuum mechanics approach, a parallel GIMP computational method is presented using the Structured Adaptive Mesh Refinement Application Infrastructure (SAMRAI). SAMRAI is used for multiprocessor distributed memory computations, as a platform for domain decomposition, and for multilevel refinement of the computational domain. Nested computational grid levels (with successive spatial and temporal refinements) are used in GIMP simulations to improve the computational accuracy and to reduce the overall computational time. The domain of each grid level is divided into multiple rectangular patches for parallel processing. This domain decomposition embedded in SAMRAI is very flexible when applied to GIMP. As an example to validate the parallel GIMP 5 computing scheme under SAMRAI parallel computing environment, numerical simulations with multiple length scales from nanometer to millimeter were conducted on a 2D nanoindentation problem. A contact algorithm in GIMP has also been developed for the treatment of contact pair between a rigid indenter and a deformable workpiece. GIMP results are compared with finite element results on indentation for validation. A GIMP nanoindentation problem with five levels of refinement was modeled using multiprocessors to demonstrate the potential capability of the parallel GIMP computation. 2.1 Introduction The material point method (MPM) has demonstrated its capabilities in addressing such problems as impact, upsetting, penetration, and contact (e.g. Sulsky, Zhou and Schreyer (1995); Sulsky and Schreyer (1996)). In MPM, two descriptions are used  one based on a collection of material points (Lagrangian) and the other based on a computational background grid (Eulerian), as proposed by Sulsky, Zhou and Schreyer (1995). A fixed structured mesh is generally used in the background throughout the MPM simulations. The material points are followed throughout the deformation of a solid to provide a Lagrangian description and the governing field equations are solved at the background grid nodes so that MPM is not subject to mesh entanglement. Compared to the finite element method (FEM), MPM takes advantage of both Eulerian and Lagrangian descriptions and possesses the capability of handling large deformations in a more natural manner so that mesh lockup problems present in FEM are avoided. Additionally, for problems involving contact, MPM provides a naturally nonslip contact algorithm to avoid the penetration between two bodies based on a common background mesh (Sulsky, Zhou and Schreyer (1995); Sulsky and Schreyer (1996)). One drawback of the 6 conventional MPM is that when the material points move across the cell boundaries during deformation, some numerical noise/errors can be generated, Bardenhagen and Kober (2004). To solve the instability problems associated with the conventional MPM simulations, Bardenhagen and Kober recently proposed the generalized interpolation material point method (GIMP) and implemented for onedimensional simulations. The present investigation extends the GIMP presented by BardenhagenKober to twodimensional simulations and applies it to simple tension and indentation problems. Furthermore, a refinement technique and a parallel processing scheme are developed so that the serial GIMP algorithm and code can be extended for parallel computation of large scale computations based on the continuum mechanics approach. Parallel processing has been used successfully in numerical analysis using different methods, such as FEM and boundary element method (BEM) (Mackerle (2003)) and molecular dynamics (MD) (Kalia and Nakano (1993)). The computational time on parallel processors can be reduced to a small fraction of the time consumed by a single processor at the same speed. Parallel processing generally involves issues, such as domain decomposition/partitioning, load balancing, parallel solver/algorithms, parallel mesh generation, and multigrid (Mackerle (2003)). Domain decomposition has been widely applied in parallel processing in FEM (Hsien (1997)). With partitioning of the overall computational domain, sequential FEM algorithm usually cannot be used directly in parallel processing without some modification, primarily due to the coupling of a large number of simultaneous linear equations. Remeshing is sometimes needed in each subdomain. The interfaces of neighboring subdomains must be meshed identically for subsequent communications (Mackerle (2003)). These problems are intrinsic to certain 7 numerical methods, such as FEM; however, they can be totally or partially avoided if other appropriate computational methods are used. For example, the domain decomposition is more straightforward for structured meshes, and large systems of coupled equations can be avoided, if explicit time integration is used. Recently a platform for parallel computation, namely, the structured adaptive mesh refinement application infrastructure (SAMRAI) (Hornung and Kohn (2002)), has been developed by the Center for Applied Scientific Computing at the Lawrence Livermore National Laboratory. SAMRAI has provided interfaces for userdefined data types so that material points carrying physical variables (mass, displacement, velocity, acceleration, stress, strain, etc.) can be readily defined. As a result, SAMRAI is very suitable for handling material points and their physical variables in MPM or its variant, GIMP. In this investigation SAMRAI is used for parallelizing GIMP. SAMRAI has also provided a foundation for parallel adaptive mesh refinement (AMR) with the use of either dynamic or static load balancing (Wissink, Hysom and Hornung (2003)). This function allows SAMRAI to process both spatial and temporal refinements in areas of interest, typically with high gradients in some physical variables (e.g., strains), and to use coarse mesh in the remaining areas. With the appropriate use of fine and coarse meshes in different regions, multiscale simulations using MPM can provide desired computational accuracy with reduced costs associated with computer memory and computational time. Material multiscale simulations span from electronic structure, atomistic scale, crystal scale, to macro/continuum scale (Horstemeyer, Baskes, Prantil, Philliber and Vonderheide (2003); Komanduri, Lu, Roy, Wang and Raff (2004)). Appropriate simulation algorithms can be used at various scales, e.g., ab initio computation for 8 electronic structure, molecular dynamics at the atomistic scale, crystal plasticity or mesoplasticity at the crystal scale, and continuum mechanics at the macro scale (Horstemeyer, Baskes, Prantil, Philliber and Vonderheide (2003)). At the continuum scale, FEM is generally used. Recently, the meshless local PetrovGalerkin (MLPG) method (Shen and Atluri (2004a, 2005)) and the continuum/lattice Green function method (Tewary and Read (2004)) have been used to couple with molecular dynamics seamlessly. The MLPG method is a simple, lesscostly alternative approach to FEM (Atluri and Shen (2002)). For the purposes of providing the insights into the discrete atomistic system and coupling with continuum, an equivalent continuum was defined in the MD region to compute the atomic stress based on the Smoothed Particle Hydrodynamics (SPH) method (Shen and Atluri (2004b)). The atomic stress tensor computed using the SPH method is more natural than other atomic stress formulations because it is in the nonvolumeaveraged form and rigorously satisfies the conservation of linear momentum. Hence, it is applicable to both homogeneous and inhomogeneous deformations. A tangent stiffness formulation was developed for both MLPG and MD regions and the displacements of the nodes and atoms are solved in one coupled set of linear equations. The MLPG/MD coupling has been demonstrated to be capable of enforcing the local balance equations in the handshaking region between continuum mechanics and molecular dynamics (Shen and Atluri (2005)). The simulation using parallel GIMP computing scheme in this investigation will focus on multiscales, e.g., from nanometer to millimeter, based on the continuum mechanics approach, namely, 2D GIMP. An example used for validating the simulation at several length scales at the continuum level is nanoindentation. It involves the contact issue 9 between a rigid indenter and a deformable workpiece. A contact algorithm, which allows the contact interface to be located in a few computational domains, is introduced in this study. The contact pressure is determined from solving a set of equations from multiple processors. Parallel GIMP results on nanoindentation are compared with FEM results using the ABAQUS/Explicit code. A nanoindentation model with five levels will be used; this model allows simulation from nanometer to millimeter scales. 2.2 Generalized Interpolation Material Point (GIMP) Method The governing equations in both conventional material point method (MPM) (Sulsky, Zhou and Sheryer (1995); Hu and Chen (2003); Bardenhagen (2002)) and generalized interpolation material point (GIMP) method (Bardenhagen and Kober (2004)) are briefly summarized in this section. The weak form of the momentum conservation equation in the conventional MPM is given by ∫ ⋅ Ω = −∫ ∇ Ω+ ∫ ⋅ + ∫ ⋅ Ω Ω Ω ∂Ω Ω ρw ad ρss : wd ρcs wdS ρw bs d , (21) where w is the test function, a is the acceleration, and ss, cs and b are the specific stress, specific traction, and specific body force, respectively. Ω is the current configuration and ∂Ω is the surface with applied traction. The material density, ρ, can be approximated as the sum of material point masses using a Dirac delta function Σ= = − Np p t p t M 1 ρ (x, ) δ (x x ) , where Np is the total number of material points and Mp is the mass of the material point. Upon discretization of Eq. (21) using the shape functions ( t ) i p N x , the governing equations at the background grid nodes become (see Sulsky, Zhou and Schreyer (1995)) ( )int ( t )ext i t i ti t i m a = f + f . (22) 10 where the lumped mass matrix is given by Σ= = N p p t p i p t i m M N 1 (x ) , (23) and the internal and external forces are given by Σ= = − ⋅∇ p t p N p i s t p p t i M N 1 ( )int , x f s , (24) Σ Σ = = = − − + p N p p t i p t p p N p t i p s t p p t ext i M h N M N 1 1 (f ) c , 1 (x ) b (x ) , (25) where h is the thickness of a boundary layer. At each time step, all variables for each material point, such as mass, velocity, and force are extrapolated to the grid nodes of the cell in which the material point resides. New nodal momenta are computed and used to update the physical variables carried by the material points. Thus, material points move relative to each other to represent deformation in a solid. A spatially fixed background grid is used throughout the MPM computation. MPM has already demonstrated its capabilities in solving a number of problems involving impact/contact/penetration. In case of large deformation, however, numerical noise, or errors have been observed, especially when material points have just crossed cell boundaries resulting in instability problems in the MPM simulations (see, e.g., Sulsky, Zhou and Schreyer (1995); Hu and Chen (2003); Bardenhagen and Kober (2004)). The primary cause for the problem has been attributed to the discontinuity of the gradient of the shape functions across the cell boundaries (see, e.g., Hu and Chen (2003); Bardenhagen and Kober (2004)). To resolve this problem, Bardenhagen and Kober (2004) proposed a generalized interpolation material point (GIMP) method. In GIMP, the interpolation between node i and material 11 point p is given by the volume averaged weighting function ∫ Ω∩Ω = p S d V S p i p ip (x) (x) x 1 χ , (26) where Vp is the current volume of the material point, (x) p χ is the characteristic function of the material point, and Si(x) is the node shape function. The role of the weighting function is the same as the shape function in conventional MPM. The modified equation of momentum conservation (Bardenhagen and Kober (2004)) can be written as Σ ∫ ∫ Σ ∫ Σ ∫ Ω∩Ω ∂Ω Ω∩Ω Ω∩Ω ⋅ + ⋅ ⋅ + ∇ = b v x c v x v x σ v x p d d V m d d V p p p p p p p p p p p p p p δ δ χ δ χ δ χ : & , (27) where δv is an admissible velocity field, p p& is the rate of change of the material point momentum. Eq. (27) can be further discretized and solved at the grid nodes, Bardenhagen and Kober (2004). Herein, the weighting function S ip is C1 continuous under the spatially fixed background grid. Consequently the noises associated with material point crossing cell boundaries in the conventional MPM can be minimized. In this chapter, the GIMP presented by BardenhagenKober is implemented for twodimensional simulations. a refinement technique and a parallel processing scheme to extend the serial GIMP algorithm to code large scale parallel computing have been developed. The capability of parallel GIMP computing has been demonstrated by modeling nanoindentation problem. A contact algorithm has been developed to address the contact problem between the rigid indenter and the deformable workpiece. Next, the contact algorithm developed in this investigation is described. 12 2.2.1 Contact algorithm in GIMP Nanoindentation involves a contact pair of a rigid indenter and a deformable workpiece. The contact interaction between these two surfaces is governed by the Newton’s third law and Coulomb’s friction law as well as the boundary compatibility condition at the contact interface (Oden and Pires (1983); Zhong (1993)). While MPM can prevent the penetration at the interface automatically, it uses a single mesh for the two bodies. At the contact surface, all components of the variables are interpolated to the nodes from both bodies using Eqs. (23)(25). As a result, MPM using a single mesh tends to induce early contact in approaching and late separation when two parts move away from each other. So, MPM cannot model the contact behavior between two parts correctly. Hu and Chen (2003) proposed a multimesh MPM algorithm to release the noslip constraint inherent in the MPM using a single mesh. In the multimesh MPM, in addition to a common mesh for all objects, there is an individual mesh for each of the objects under consideration. All meshes are identical, i.e. nodal locations are the same. The multimesh can be generated by creating multiple nodal fields for each node. Each nodal field corresponds to an object. In multimesh MPM scheme, the nodal masses and forces are mapped from the material points of each object to its own mesh. The nodal values are transferred to the corresponding nodes in the common mesh. When the values at a node of the common background mesh involve contributions from two parts, the contact between two parts occurs so that this node is defined as an overlapped node. Otherwise, two parts move independently. This multimesh algorithm can handle sliding and separation for the contact pair. However, in using the multimesh for contact problem in GIMP, the interaction at the overlapped nodes is still activated too early before the actual contact of 13 the material points occurs. Fig. 21 is an example illustrating early contact when Part 1 is moving toward Part 2. The four bottom particles of Part 1, labeled in hollow circles, have come into the cells of Part 2, the nodes of which cells are labeled in three dashed circles. Physical variables (e.g., normal force, and velocity component normal to the contact surface) in Part 1 will be interpolated onto these three overlapped nodes. The physical variables in the three overlapped nodes will be further interpolated into material points within the top layer of cells in Part 2, and contribute to the stress and deformation in the entire Part 2. With this treatment in the previous multimesh algorithm, even though Parts 1 and 2 are not in physical contact, the particles in Part 2 will contribute to the physical variables of particles in Part 1, leading to numerical early contact, and vice versa, through the overlapped nodes. Similar situation occurs when Part 1 is retracting from Part 2, resulting in late separation of two parts. Unless other measures are taken to prevent these physically incorrect early contact and late separation problems, they could cause large errors in GIMP and must be corrected in contact problems. Fig. 21: Illustration of early contact in multimesh In this chapter, a new contact algorithm is developed for GIMP simulations. Fig. 22 illustrates the contact algorithm for the contact pair between a rigid indenter and a Part 1 Part 2 Overlapped nodes 14 deformable workpiece. Although circular points are used in this schematic diagram, it should be noted that the points are representations of areas occupied by these points, based on the GIMP algorithm. A frictionless contact is assumed in this investigation. At the beginning of a time step, a material point is located at point A. At the end of this time step, the material point moves to B, if there is no contact interaction. Fig. 22: Schematic of contact algorithm between the rigid indenter and the deformable workpiece To satisfy the displacement compatibility condition, the material point has to be brought to the indenter surface and kept in contact with the indenter. The contact velocity correction c p V can be determined based on the rigid surface orientation indicated by its unit outward normal vector n. The final location of the material point is set to C by a contact pressure. Hence, the velocity of a material point p under contact can be determined by Σ Σ Σ Σ Δ + + Δ = + + Δ = = N i ip i c i N i ip i i i N i ip i c i i i N i p i ip S m S t m t S m S t p F F V V p F F 0 0 0 ( 0 ) , (28) where N is the number of nodes contributing to material point p, c i F is the contact force on node i, 0 i p and 0 i F are the nodal momentum and force without consideration of the A C B n D 0 p V c p V p V Rigid surface Rigid indenter Indenter side Workpiece side Workpiece 15 contact, respectively. The velocity 0 p V of the material point without the consideration of contact is given by Σ + Δ = N i ip i i i p S m 0 0 t V0 p F . (29) The contact force c i F on node i is the resultant of the contact pressures on the neighboring particles, and can be computed in terms of contact pressures using the approach given by Bardenhagen and Kober (2004), i.e., Σ ∫ = ∂Ω = Q q c i q c i q S ds 1 F (x)P , (210) where Q is the total number of material points in contact with the indenter. If the contact pressure c q P is assumed to be constant in the contact area occupied by material point q, then, Σ= = Q q c iq q c i T 1 F P , where ∫ ∂Ω = q T S ds iq i (x) . Since c p p p V = V0 + V , the contact velocity c p V for material point p is given by = Δ Σ Σ N i Q q c iq q i c ip p T m V t S P . (211) Eq. (211) can be established for each material point in contact. At each material point there is an unknown contact pressure c q P . Therefore, the number of unknown pressures, c q P , is equal to the number of points in contact. In parallel computing, points in contact might be located in different domains processed by different processors. Consequently, a parallel solver is needed to solve Eq. (211) in this investigation. Since contact can only occur on the outer surface of an object, Eq. (211) is solved analytically under the physical contact condition Pc ⋅n < 0 q to find the contact pressure at all material points in 16 contact with the indenter. The contact pressure is then extrapolated to the nodes from the contact material points to update the total nodal forces. 2.2.2 Numerical implementation Fig. 23: Material points in cells and the weighting function in 2D GIMP Considering the case where initially there are four material points in a cell, the 2D weighting function is depicted in Fig. 23. To compute the weighting function, χp( x ) is one in the current region occupied by the material point p and zero elsewhere. In this figure, one node is at the origin and the horizontal axes give material point positions normalized by the cell size. Fig. 23 is based on the same material point characteristic function and node shape function as in Bardenhagen and Kober (2004). It is noted that the computation of the weighting function in the deformed state involves some practical difficulties because the integration boundaries in Eq. (27) can be difficult to obtain. To circumvent this problem, it is assumed that the shape of the region occupied by the four material points remains rectangular without rotation, so that Eq. (26) can be evaluated analytically. This assumption leads to significant saving in the computational time while introducing only small errors. Using this assumption, GIMP is extended to 2D simulations and the results are presented in Section 4. 1 0.5 0 0.5 1 1.5 1.5 1 0.5 0 0.5 1 1.5 0 0.2 0.4 0.6 0.8 17 2.3 Parallel Computing Scheme Using GIMP with SAMRAI 2.3.1 SAMRAI The Structured Adaptive Mesh Refinement Application Infrastructure (SAMRAI), a scientific computational package for structured adaptive mesh refinement and parallel computation, is used with the GIMP for parallel computation of largescale simulations. SAMRAI is chosen because of its similarity in grid structure with GIMP. In GIMP, the computation is usually independent of the background grid mesh so that a structured spatiallyfixed mesh can be used throughout the entire simulation process. This advantage makes GIMP highly suitable for parallel computation, as the domain decomposition for structured mesh can be easily performed and no remeshing is required. Thus the complexity and inefficiency associated with parallel processing can be avoided. In SAMRAI, the computational domain is defined as a hierarchy of nested grid levels of mesh refinement (Berger and Oliger (1984)) as shown in Fig. 24. Each grid level is divided into nonoverlapping, logicallyrectangular patches, each of which is a cluster of computational cells. Indices are used extensively in SAMRAI to manage grid levels and patches. For example, patch connectivity is managed by the cell indices. The organization of the computational mesh into a hierarchy of levels of patches allows data communication and computation to be expressed in geometricallyintuitive box calculus operations. Communication patterns for data dependencies among patches can be computed in parallel without interprocessor communications, since the mesh configuration is replicated readily across processor memories. Interprocessor communications, i.e., data communications between patches on the same as well as neighboring levels, are predefined by SAMRAI communication schedules. Problem 18 specific communication interfaces are also provided by SAMRAI. SAMRAI supports several data types defined in a patch, such as cellcentered data, nodecentered data, and facecentered data. These data are stored as arrays to allow numerical subroutines to be separated easily from the implementation of mesh data structures. Userdefined data structures over a patch, which can be accessed through cell index, are supported by SAMRAI. These characteristics make SAMRAI a very flexible parallel computing environment for numerous physics applications (Wissink, Hysom, and Hornung (2003)). Fig. 24: Illustration of a hierarchy of three nested grid levels of mesh refinement 2.3.2 Spatial and temporal refinements In the application of SAMRAI to largescale GIMP simulations, the techniques for refinement, both spatial and temporal, have to be developed to achieve high accuracy in areas of high stress/strain gradients while reducing the overall computational time by using coarse mesh in regions of low stress/strain gradients. Since a structured mesh is used in GIMP, the refinement can be implemented by imposing fine levels of subgrids at locations of interest, using the approach adopted by Berger and Oliger (1984) in SAMRAI. The scheme for the structured grid refinement is illustrated in Fig. 24. The Level 1 Level 2 Level 3 19 cell size ratio, also called the refinement ratio, of two neighboring levels is always an integer for convenience. The advantage of this refinement technique is that nesting relationships between different levels can be handled. A material point in GIMP can be split into several small material points. Tan and Nairn (2004) proposed a criterion to split material points based on local deformation gradient. If the refinement ratio is two in each direction, one coarse material point can be split into four material points in the next fine level in 2D case. However, this splitting technique can become complicated when conservation of energy and momentum have to be enforced. In this chapter, a more natural refinement approach is developed to avoid direct splitting and merging processes by using material points of the same size and mass in overlapped region (called ghost region) between two neighboring levels. Fig. 25 shows two neighboring coarse and fine grid levels in 2D GIMP computations with a refinement ratio of two. The thick line represents the physical boundary of the fine level with four layers of ghost cells. Initially, four material points are assigned to each cell at the fine level. At the coarse level, the portion overlapped by the fine level is assigned with 16 material points per cell. Hence, these material points have the same size and initial positions as those at the fine level. The rest of the coarse level is assigned with four material points per cell. GIMP provides a natural coupling of the material points with different sizes at the same grid level. This is because the weighting function depends on the characteristic size of the material points and cell length and the interpolation between the nodes and material points is weighed by the mass of the material point. In GIMP computation, each level is computed independently with the physical variables communicated through the ghost regions between neighboring levels. Two data exchange 20 processes, namely, refinement and coarsening are used in the communication. Refinement process passes information from the coarse level to the immediate fine level, while coarsening process will pass information from fine level to the next coarse level. In the refinement process, physical variables at the fine material points inside the thick lines in Fig. 25 are copied directly to replace the material points at the coarse level. In the coarsening process, the physical variables at coarse material points are copied to the ghost cells of the immediate fine level. In the refinement, the material points located in the ghost cells at the fine level are eliminated first, and the material points in the corresponding region of the coarse level (with the same size as points in the immediate fine level) are copied to ghost cells at the fine level. In the copying process, if some (small) material points in the coarse level fall into the interior (inside the square on the right of Fig. 25) of the immediate fine level, these points will not be copied to this fine level, as points in the fine level will be able to carry over all computations in the interior of the fine level already. In coarsening, the material points of the coarse level located in the region overlapped by the fine level (inside the square on the left of Fig. 25) are eliminated first, and the material points in the fine level are then copied to the immediate coarse level. With these refinement/coarsening operations, the material points can move around freely, including moving outside the original level in large deformations. To ensure that this coarsening process can still be performed reliably during deformation, sufficiently wide region of cells should be assigned with refined material points at the coarse level so that the ghost cells of the fine level always stay within the region with fine material points on the coarse level. At the coarse level, the interior cells covered by the fine level do not participate in 21 the computation and there are no material points inside (see Fig. 25). Fig. 25: Two neighboring coarse and fine grid levels in 2D GIMP computations Fig. 26: Nested multilevel refinement and reduction in the number of cells with number of levels The refinement techniques can be applied for multiple times at the regions of interest, such as the stress concentration regions. A fixed refinement ratio of two between two neighboring levels is very effective in reducing the total number of computational cells. Fig. 5 shows nested multilevel refinement and its corresponding relation between the total number of cells and the number of grid levels. The cell percentage represents the ratio of the total number of cells with multilevel refinement mesh to the total number of cells with onelevel finest mesh. If each fine level occupies one quarter of the neighboring coarser level, as shown in Fig. 26 (a), the cell percentage as a function of Coarse Fine Level 1 Level 2 3 (b) 0 20 40 60 80 100 1 2 3 4 5 6 7 Number of levels Cell percentage (%) (a) 22 the number of grid levels can be calculated, as shown in Fig. 26 (b). For example, when totally four levels of successive refinements are used the total number of cells is about 8% of that of one uniform fine mesh. A reduction in the number of computational cells leads to a reduction in the number of material points. Hence, the total amount of computational time can be reduced significantly. However, refinement and coarsening communications will cost additional computational time, as will be discussed in Section 4. Another advantage of the multilevel refinement is that it allows for temporal refinement. Since the computation at each grid level is conducted independently, different time step increments can be used for computation at different levels. For example, a smaller time step increment can be used for the fine level to improve computational accuracy, while a larger time step increment can be used for the coarse level. Since the refinement ratio is an integer, the time step increment ratio should also be an integer for convenience in the computation and data communication/synchronization. For example, in Fig. 25, when the refinement ratios in both directions are fixed at two, the time step increment ratio should be set to two as well. As a result, two time step computations are performed at the fine level, and results are passed over to the immediate coarse level to couple with the results at the coarse level. 2.3.3 Domain decomposition GIMP uses structured mesh, consistent with SAMRAI, so that domain decomposition is straightforward and no remeshing, in general, is necessary. Fig. 27 (a) shows a twodimensional computational domain decomposed into two patches separated by a horizontal dash line. The elliptical solid object with different boundary conditions applied at different regions is inside this domain/grid. After discretization, there are a certain 23 number of material points and part of the boundary in a patch, which will be computed individually. It may be noted that patch boundary does not have to coincide with the boundary of the material continuum. The patch boundary is always chosen to be larger than the region occupied by the material continuum so that there is extra space for the material to deform. This will not cause any additional computational burden as the GIMP computation is only carried out on material points inside the patch. Each patch can be processed by a single processor and the convenience in creating patches will provide great flexibility in parallel processing. Fig. 27: A computational domain of two patches in one grid level Communication between two neighboring patches is realized through information sharing in the region overlapped by the two patches. The overlapped regions are also called “ghost” regions, as shown in Fig. 27 (b). The ghost cells are denoted by dash lines. For ease of visualization, only the ghost cells overlapped by the other patch are shown and the ghost cells along the other three sides of a patch are not shown. On one grid level, patches can communicate with each other by simply copying data from one patch to another at the same computation time (Fig. 27 (b)). Using the material point information from the previous time step, and the physical boundary conditions, each patch is ready to advance one more time step. At this time, the material point information in the outermost Patch 2 Patch 1 Ghost of patch 1 Copy from patch 2 Ghost of patch 2 Copy from patch 1 (a) (b) 24 layer of the ghost cells becomes inaccurate. For instance, one outermost grid node in patch one, marked by the circle, obtains information from eight material points before advancing to the next step. After advancing, it extrapolates to eight material points. However, in patch two, the grid node at the same location obtains information from sixteen surrounding material points. It extrapolates to these sixteen material points after advancing. Typically, after each step, the material points in the next inner layer in the ghost region become inaccurate as well. Ghost cells and material points are attached to each patch to ensure accuracy of the interior. Each patch can be computed independently for one GIMP step since the momentum conservation equation is solved at each node and there are no coupled equations to solve. No data exchange is necessary during the GIMP step. Therefore, different patches can be assigned to different processors for parallel processing. After one GIMP step, the data in the ghost cells will be updated. Fig. 28: Flowchart showing advancement of grid levels recursively starting from the coarsest to the finest level in GIMP AdvanceLevel (level_number) Advance each patch in level by level time step Start level_number=0 Next finer level exists AND level time > next finer level time? level_number = level_number+1 Data communication for level YES NO End 25 Copying material points to ghost cells involves data exchange between processors, which costs additional time. The more the number layers of ghost cells, the longer the time needed for communication, but communication can be performed less frequently. A minimum of two layers of ghost cells are necessary to ensure that computation at the material points inside a patch is always correct. If three levels of ghost cells are chosen, the communication can be performed after every two increments of each patch. With these refinements and domain decomposition schemes for GIMP, it is possible to implement GIMP into the SAMRAI platform. In this study, the refinement ratio is chosen as two. Four layers of ghost cells are augmented to a patch such that data communications, including both data exchange on the same level and between neighboring levels, are performed every two timesteps for each fine grid level. This is critical because data exchange between levels has to be performed when the two levels are synchronized. Fig. 28 shows the flowchart advancing all grid levels recursively starting from the coarsest level for one coarsest time step. It may be noted that the sequential GIMP algorithm can be used to advance each patch without modification. 2.4 Numerical Examples and Results A Beowulf Linux cluster of 8 identical PCs were used in the simulations. Each PC has a Pentium 4 processor with a 2.4 GHz CPU, 512 MB RAM except that the master node has a memory of 1 GB. A gigabit switch is used to connect the network. Two examples are used for validation of the 2D parallel GIMP computing under SAMRAI platform. The first example is simple tension of polycrystalline silicon under plane strain conditions. The material is assumed to be homogeneous, isotropic, and linear elastic. The Young’s modulus is 170 GPa and the Poisson’s ratio is 0.18. One end is 26 constrained along the X direction while a normal traction is ramped up on the other end. The size of the tensile model is 0.06 mm × 0.04 mm. The length of a square grid cell is 0.002 mm and the time step is 5×104 ps. For verification, the same problem was simulated using both conventional MPM and FEM (ABAQUS/Explicit). Fig. 29 shows GIMP, MPM and FEM simulation results of normal stresses in the Xdirection at different increments from a simple tension problem. The simulation using the conventional MPM in Fig. 29 (a) shows material separation close to the free end with severe numerical instability after 275 increments. Fig. 29 (b) and (c) show the normal stress distribution in the tensile direction and deformation after 500 increments from FEM and GIMP. It may be noted that FEM results show a stress contour plot on the deformed mesh while the GIMP results show a discrete scattered plot of material points. These two results are in good agreement with the difference in the maximum value being less than 10%. Fig. 29: Simulation results of tensile stress contours for a simple tension problem (c) FEM (a) Conventional MPM (b) GIMP 27 Fig. 210: Loading conditions for a simple indentation problem In the second example, indentation on the same silicon material is simulated. The workpiece is subjected to a pressure applied in the middle of the top surface (Fig. 210 (a)) under plain strain conditions with a thickness of 0.001 mm for computing the mass and forces. The magnitude of the pressure increases linearly with time for the first 1500 increments, and is then kept constant (see Fig. 210 (b)). The cell size is 0.001 mm in both directions and the time step is 20 ps for both FEM and GIMP simulations. Due to symmetry, only one half of the workpiece is modeled. This simulation is performed with two patches in one uniform grid level. Two processors are used and one patch is assigned to each processor. Fig. 211 shows GIMP and FEM results of normal stresses in the Ydirection at different increments. The dashed line in Fig. 211 (a) is the boundary between the two patches. Fig. 211 (a) and (b) are plots of normal stresses in Ydirection at 500 time increments for GIMP and FEM simulations. The difference in stress values in Fig. 211 (a) and (b) is less than 5%. It should be noted that the FEM simulation aborted at 1348 increments due to excessive element distortion. The GIMP simulation did not encounter this problem. Fig. 211 (c) shows the GIMP stress result after 2000 increments. This demonstrates the capability of GIMP in handling excessive distortions. X Y (a) (b) Time Pressure t = 3×10−8 s 60 GPa 28 Fig. 211: GIMP and FEM results of normal stress variation in the Ydirection at different increments In order to validate the multilevel refinement algorithm and parallel communication, as well as the proposed contact algorithm, a simulation of nanoindentation with a wedge indenter was conducted under 2D plane strain conditions. The workpiece is aluminum and the indenter is assumed to be rigid. Fig. 11 shows the indentation model. The area below the indenter where high stress gradients are expected is refined, as shown in Fig. 212 (a). A prescribed velocity was applied on the indenter, as shown in Fig. 212 (b). (c) GIMP at 2000 increments (a) GIMP at 500 increments (b) FEM at 500 increments 29 The work piece dimensions are 60 μm × 40 μm. It is fixed in the Ydirection at the bottom. Only half of the model is simulated because of symmetry. The cell sizes are 500 nm, 250 nm and 125 nm for levels 1, 2 and 3, respectively. Each level is divided into four patches with approximately the same size. The maximum indentation depth in the simulation is about 450 nm. The dotted lines in Fig. 212 (a) illustrate the four patches in level 1. For comparison, an explicit FEM simulation (using ABAQUS/Explicit) was carried out under the same conditions. The FEM element size is uniform and is the same size as the finest GIMP background grid size. In this example, the maximum indentation depth (450 nm) is relatively small compared to the finest cell size, so that FEM simulation has not encountered excessive mesh distortion. Fig. 212: Schematic of 2D indentation showing (a) three levels of refinement and (b) the indenter velocity history Fig. 213 shows a comparison of contours of normal stresses in the Ydirection at the maximum depth for both FEM and parallel GIMP simulations. The axis of symmetry of the workpiece is located at X = 0.03 mm. For FEM, the plot is the contour of nodal stresses with deformed positions, and for GIMP, it is a discrete scattered plot of stress at deformed material points. The area below the indenter with high stress gradients is refined as shown in Fig. 212 (a) for parallel GIMP computation. The borders of grid levels 2 and 3 can barely be seen in Fig. 213 (b) due to the use of high material point (a) (b) Time V 20 m/s V Level 1 Level 2 3 X Y 30 density. Fig. 214 is a closeup view of shear stresses in which the three grid levels are shown. Results show that the normal and shear stresses from both parallel GIMP and FEM simulations agree very well. The difference of the normal stresses in the Ydirection for the material point in contact with the indenter tip and the stress of the FEM node at the same location is 4.4%. It may be noted that some nonsmoothness in the GIMP stresses around the level boarders can be seen. This nonsmoothness is caused by the refinement and coarsening and the error associated with this is negligible for these simulations. Fig. 213: Comparison of normal stresses in Ydirection from FEM and GIMP GIMP simulations using a uniform cell size of 500 nm and 125 nm were performed under the same conditions as in Fig. 212 to further verify the refinement/coarsening algorithm. Fig. 215 shows normal stresses in the Ydirection and shear stresses in GIMP simulations using 500 nm uniform cells. In this case the material points in contact with the indenter are 16 times (4 times in each direction) larger than those with two levels of refinements. In general, the stress magnitudes agree with those in Fig. 213 and Fig. 214. Fig. 216 (a) shows indentation load versus depth curves from FEM and GIMP with (a) FEM (b) GIMP 31 different grid sizes. From GIMP computations, the load versus depth curves with three levels of refinement agree very well with the results from a uniform finest mesh. The load versus depth curve from the FEM simulation with a uniform cell size of 125 nm under the same boundary conditions is plotted for comparison. It can be seen from Fig. 216 (a) that the trend of the load versus depth plots from FEM and GIMP simulations are similar. The difference in indentation load at the end of loading, which corresponds to 450 nm of indentation depth, is 5.9% between FEM and GIMP with 3 grid levels. When the depth is less than 100 nm, there is only one material point in contact with the rigid indenter. The assumption of constant pressure causes large differences under this circumstance. However, if the GIMP cell size is further refined to 62.5 nm and the size of the material point is 31.25 nm, the difference between GIMP and FEM becomes smaller, as can be seen in Fig. 216 (b). Fig. 214: Comparison of shear stresses from FEM and GIMP Other simulations were conducted for the same problem with three levels of refinement using different number of processors to test the efficiency of parallel computing. The number of patches at each level is the same as the number of processors and the size of (a) FEM (b) GIMP 32 each patch is approximately the same. The resultant stress distribution and indentation load versus depth plots are the same as the previous results. The average time per computational step is 7.14 s when one processor is used and is reduced to 4.26, 3.40, 2.18 s, respectively when two, three, and four processors are used. When four processors are used, the CPU time per step is only 30.5% of that of one processor. This gives a speedup by a factor of 3.28. In the ideal case without communication overhead, the speedup would be 4. The reduction in speedup from the ideal number is because of the time involved in data communication between processors. It has been observed that the refinement and coarsening algorithm consume most of the communication time. Moreover, in refinement and coarsening, most of the time is taken to search for the corresponding material point in another grid level. This portion of the computational time can be reduced, if improved searching algorithm or more optimized algorithm for the storage of material points can be implemented. Fig. 215: Normal and shear stresses of GIMP simulations with a uniform cell size of 500 nm The manual refinement for the indentation problem is adequate since the region of high (a) Normal Stress (b) Shear Stress 33 stress gradient is known to occur below the indenter. The finest level covers the indenter and part of the specimen. With the same initial condition, the results at the finest level is identical to the results in the same area if a uniform fine mesh is used for the entire domain that requires much longer computational time. The computational load of each processor is balanced statically by assigning approximately the same number of material points to each processor. Dynamic load balance is supported by SAMRAI and can potentially improve the efficiency of the simulation. Fig. 216: Indentation load versus depth curves from FEM and GIMP with different grid sizes To demonstrate the capability of the algorithm developed in this investigation for multiscale simulation, an indentation model with multiple length scales is simulated with eight processors. The dimensions of the workpiece are 0.25 mm × 0.125 mm. Initially, the velocity of the indenter increases from 0 to 150 m/s linearly with time and is then kept constant. Five successive levels of refinement are used in this simulation. The smallest material point represents an area of 64 nm × 64 nm, and the largest material point covers an area of 1 μm × 1 μm. Each level is divided into 8 equalsized patches for best load balance. Since the contact surface can evolve into several patches, a parallel 0 1 2 3 4 5 6 0 100 200 300 400 500 Depth (nm) Load (mN) GIMPMcoarse (500 nm) GIMPMfine (125 nm) GIMPM3 levels FEMfine (125 nm) 0 0.5 1 1.5 2 2.5 0 40 80 120 160 Depth (nm) Load (mN) GIMPM500 nm GIMPM125 nm GIMPM62.5 nm FEM62.5 nm (a) (b) 34 solver is implemented to solve Eq. (211) to find the contact pressure based on the Portable, Extensible Toolkit for Scientific Computation (PETSc). An aluminum workpiece is chosen with the Young’s modulus and Poisson’s ratio of 70 GPa and 0.33, respectively. The maximum indentation depth was 9.8 μm in this simulation (i.e., 153 times the size of the finest material point). It took nine hours to simulate this problem with eight processors. Fig. 217 (a) gives the normal stress distributions and Fig. 217 (b) shows normal stress distribution for the finest two levels. The relative large deformation in this multiscale nanoindentation problem could not be handled by FEM due to excessive distortion in the FEM mesh. However, the parallel GIMP code was able to complete the entire loading/unloading processes without any difficulty. This example shows clearly the advantage of GIMP for multiscale simulations over FEM. Fig. 217: Multiscale simulation of nanoindentation with five levels of refinement 2.5 Conclusions The following are specific conclusions based on the results of this investigation: A 2D generalized interpolation material point (GIMP) method has been implemented to address problems, such as particle flyingoff and alternating stress sign associated with (a) (b) 35 conventional MPM in case of relatively large deformation. To conduct multiple length scale simulations, a parallel computing scheme has been presented using GIMP under SAMRAI parallel computing environment in which multilevel grids are used for spatial and temporal refinements. A refinement/coarsening algorithm, based on material points of GIMP in two grid levels, has been developed for communication between neighboring grid levels of different refinements. With increase in the refinement levels, as well as decrease in the time step increments, the computational accuracy is greatly improved in the region of interest while the overall computational time is reduced. The computation at each grid level is performed recursively to ensure that the refinement and coarsening are performed when the two neighboring levels are synchronized. 2D MPM and GIMP were applied to simple tension and indentation problems to validate the GIMP algorithm. GIMP results agree very well with FEM results for these two examples provided that the deformations are small. The noise and instability problems present in conventional MPM are not observed in the GIMP simulations. As the deformation is increased, GIMP continued to execute while FEM aborted due to element distortion. Also GIMP results are stable. Thus GIMP is able to handle relatively large deformation problems. For the nanoindentation problem, a GIMP algorithm for the contact between a rigid indenter and a deformable workpiece was developed. A reasonably good agreement between GIMP and FEM results was reached, validating the contact algorithm presented in this investigation. Another nanoindentation example with multiple length scales from a few nanometers to 36 submillimeters was simulated and numerical results validated the parallel GIMP computing with the use of SAMRAI. 37 Chapter 3 Structured Mesh Refinement in GIMP Method for Simulation of Dynamic Problems The generalized interpolation material point (GIMP) method, recently developed using a C1 continuous weighting function, has solved the numerical noise problem associated with material points just crossing the cell borders, so that it is suitable for simulation of relatively large deformation problems. However, this method typically uses a uniform mesh in computation when one level of material points is used, thus limiting its effectiveness in dealing with structures involving areas of high stress gradients. In this chapter, a spatial refinement scheme of the structured grid for GIMP is presented for simulations with highly localized stress gradients. A uniform structured background grid is used in each refinement zone for interpolation in GIMP for ease of generating and duplicating structured grid in parallel processing. The concept of influence zone for the background node and transitional node is introduced for the mesh size transition. The grid shape function for the transitional node is modified accordingly, whereas the computation of the weighting function in GIMP remains the same. Two other issues are also addressed to improve the GIMP method. The displacement boundary conditions are introduced into the discretization of the momentum conservation equation in GIMP, and a method is implemented to track the deformation of the material particles by tracking the position of the particle corners to resolve the problem of artificial separation of material particles in GIMP simulations. Numerical 38 simulations of several problems, such as tension, indentation, stress concentration and stress distribution near a crack (mode I crack problem) are presented to validate this refinement scheme. 3.1 Introduction The material point method (MPM) uses a collection of material points, mathematically represented by Dirac delta functions to represent a material continuum (Sulsky, Zhou, and Schreyer (1995); Hu and Chen (2003); Guilkey and Weiss (2003)). A spatially fixed background grid, and interpolation between grid nodes and material points are introduced to track physical variables carried by the material points in the Lagrangian description. Field equations are solved on the background grid in the Eulerian description. Physical variables are interpolated from the solutions on the background grids to material points back and forth for solution and convection of physical variables. In general, the isoparametric shape functions, same as those used in the finite element method (FEM), are used. As the MPM simulation is independent of the background grid, a structured grid is usually employed for purposes of simplicity. The movement of the material points represents the deformation of the continuum. MPM has been demonstrated to be capable of handling large deformations in a natural way (Sulsky, Zhou, and Schreyer (1995)). However, primarily due to the discontinuity of the gradient of the interpolation function at the borders of the neighboring cells, artificial noise can be introduced when the material points move just across the grid cell boundaries, leading to simulation instability for MPM. The generalized interpolation material point (GIMP) method, introduced by Bardenhagen and Kober (2004) can resolve this problem. In GIMP a C1 continuous interpolation function is used and each material point/particle occupies a region. GIMP 39 has been demonstrated to be stable and capable of handling relatively large deformations (Ma, Lu, Wang, Roy, Hornung, Wissink, and Komanduri (2005)). The current MPM typically uses a uniform background mesh for solving the field equations. However, this is not efficient when stress gradients are high such as stress concentrations in a plate with a hole, or the stress field of a workpiece under indentation. In contrast, transitional mesh is effective in solving problems involving rapidly varying stress in an area. Wang, Karuppiah, Lu, Roy, and Komanduri (2005) have presented a method using an irregular background mesh to deal with problems involving rapidly varying stress, such as stress field near a crack. However, this approach does not use regular structured background mesh so that mesh generation encounters the same difficulty as FEM, and leads to the loss of some advantage of MPM on the ease of generating mesh for a complex problem. The use of structured grid in GIMP has facilitated the implementation of GIMP in parallel processing. A refinement scheme based on splitting and merging material particles was proposed by Tan and Nairn (2002). Recently, a multilevel refinement algorithm has been developed for parallel processing using the structured adaptive mesh refinement application infrastructure (SAMRAI) (Hornung and Kohn (2002); Ma, Lu, Wang, Roy, Hornung, Wissink and Komanduri (2005)). The computational domain is divided into multiple nested levels of refinement. Each grid level is uniform but has a different cell size. Smaller material particles and smaller cell sizes are used in each finer level. Two neighboring levels are connected by overlapped material particles of the same size and data communication between levels is performed at predefined intervals. However, the refinement through material particles requires extra communication and 40 simulation time. In this chapter, a refinement for GIMP based on the transitional grid nodes is developed. This refinement is natural and does not involve extra simulation time. Moreover, the refined grid remains uniformly structured in each refinement region. While the problem associated with artificial noise has been resolved with the use of GIMP method, it has been observed recently that material separation could occur if the deformation of the material particles was not tracked (Guilkey (2005)). Tracking the deformation of material particles properly in GIMP is necessary especially when the material particles are stretched. In this chapter, an approach is developed for tracking the particle deformation to resolve material point separation problem. This chapter focuses on the refinement scheme for structured grid. Several numerical problems, such as tension, indentation, stress concentration and stress distribution near a crack (mode I crack problem) were simulated to verify this refinement algorithm, as well as to demonstrate the effectiveness of tracking particle deformations. 3.2 GIMP For the purpose of completeness, the basic equations in GIMP (Bardenhagen and Kober (2004)) are summarized here. In dynamic simulations, the mass and momentum conservation equations are given by + ρ∇⋅ v = 0 ρ dt d , and (31) ρa = ∇⋅σ + b in Ω, (32) where ρ is the material density, a is the acceleration, σ and b are the Cauchy stress and body force density, respectively. The displacement and traction boundary conditions are given as 41 u = u on u ∂Ω , (33) τ = τ on τ ∂Ω , (34) where ∂Ω ⊂ ∂Ω, u ∂Ω ⊂ ∂Ω τ and ∂Ω ∩∂Ω = 0. u τ In variational form, the momentum conservation equation can be written as ∫ ∫ ∫ ∫ Ω Ω Ω ∂Ω ⋅ = ∇ ⋅ + ⋅ − − ⋅ u ρa δvdx σ δvdx b δvdx α (u u) δvdx , (35) whereδv is an admissible velocity field, α is a penalty parameter introduced herein to impose the essential boundary conditions and α >>1 (Atluri and Zhu (1998); Atluri and Zhu (2000)). Applying the chain rule, (∇⋅σ) ⋅δv = ∇⋅ (σδv) − σ :∇δv , and the divergence theorem, Eq. (35) can be written as ∫ ∫ ∫ ∫ ∫ ∫ ∂Ω ∂Ω Ω Ω Ω ∂Ω + ⋅ − − ⋅ ⋅ + ∇ = ⋅ + ⋅ u u dS dS d d d dS u τ v u u v a v x σ v x b v x τ v δ α δ ρ δ δ δ δ τ ( ) : , (36) where u τ is the resultant traction due to the displacement boundary condition on u ∂Ω . In GIMP, the domain Ω is discretized into a collection of material particles, with p Ω as the domain of particle p. The physical quantities, such as the mass, stress and momentum can be defined for each particle. For example, the momentum for particle p can be expressed as ∫ Ω = p d p p p ρ (x)v(x)χ (x) x , where v( x ) is the velocity and ( ) p χ x is the particle characteristic function. The momentum conservation equation can be discretized as Σ ∫ Σ ∫ Σ ∫ Σ ∫ Σ ∫ Σ ∫ ∂Ω ∩Ω ∂Ω ∩Ω ∂Ω ∩Ω Ω∩Ω Ω∩Ω Ω∩Ω + ⋅ + ⋅ − ⋅ ⋅ + ∇ = ⋅ p p u p p p p p p p p p p p p p u p u p p p p d dS dS d V m d d V τ v x τ v u  u v v x σ v x b v x p δ δ α δ δ χ δ χ δ χ τ ( ) : & . (37) 42 where ∫ Ω∩Ω = p V d p p χ (x) x is the particle volume. Introducing a background grid and the grid shape function (x) i S that satisfies partition of unity 1 ) ( = Σi i S x , the admissible velocity field can be represented by the grid nodal data as =Σ i i i δv δv S (x) . Without the loss of generality, take u in Eq. (37) to be the displacement of the boundary particles at the current time step and u p u τ = σ n where u n is the unit outward normal to u ∂Ω . The momentum conservation, Eq. (37), can eventually be written for each node i as u i i b i i i p& = f int + f + f τ + f , (38) where the time rate change of nodal momentum =Σ / Δ , p i ip p p& S p t the nodal internal force vector int = −Σ ⋅∇ , p i p ip p f σ S V the nodal body force vector =Σ p p ip b i f m bS and the nodal traction force vector Σ ∫ ∂Ω ∩Ω = p i i p S dS τ f τ τ (x) . u i f is the force vector induced by the essential boundary condition given by Σ ∫ Σ ∫ ∂Ω ∩Ω ∂Ω ∩Ω = − p p i p p u i u i u p u p f σ n S (x)dS α (u  u)S (x)dS . (39) Sip is weighting function between particle p and node i given as ∫ Ω∩Ω = p S d V S p i p ip (x) (x) x 1 χ . (310) The weighting function in GIMP is C1 continuous and satisfies partition of unity. The momentum conservation, Eq. (38), can be solved at each node to update the nodal momentum, acceleration, and velocity. These updated nodal quantities can be 43 interpolated to the material particles to update the particles, as given by Bardenhagen and Kober (2004). It may be noted that the mass of each material particle does not change, so that the mass conservation equation is satisfied automatically. In the discretization of the weak form of the momentum conservation equation, a background grid is used. However, the computation is independent of the grid from one increment to another. Hence a spatially fixed structured grid can be used for convenience. In the background grid, no nodal connectivity is required and the integration is never performed on the element domain. Similar characteristics have been reported for other meshless methods, such as the meshless local PetrovGalerkin (MLPG) method (Atluri and Shen (2002)). For a uniform structured grid, the grid shape function in 3D is defined as the product of three nodal tent functions (Bardenhagen and Kober (2004)) S ( ) S (x)S ( y)S z (z) i y i x i i x = , (311) in which the nodal tent functions are in the same form, e.g., ⎪ ⎪ ⎩ ⎪ ⎪ ⎨ ⎧ ≤ − − − ≤ − ≤ + − − ≤ − ≤ − ≤ − = x i i x i x i x x i i x x i L x x x x L x x L x x L L x x x x L S x 0 1 ( ) / 0 1 ( ) / 0 0 ( ) . (312) Fig. 31 shows one 2D grid cell with four nodes. In this chapter, the particle characteristic function of the material particle located at (xp, yp) is taken as ( ) (x) y ( y) p x p p χ x = χ χ , (313) where ( x ) H[ x ( x l )] H[ x ( x l )] p x p x xp χ = − − − − + and H denotes the step function. 44 Fig. 31: 2D representation of a particle and a grid cell 3.3 Structured Mesh Refinement In this section, a refinement scheme for a structured mesh in GIMP is described. Since GIMP shares some characteristics with meshless methods, the GIMP method is expected to have h convergence in simulation (Atluri and Shen (2002)). The momentum conservation equation is essentially solved at each node (see Eq. (38)). Therefore, the number of equations to be solved is the same as the number of nodes. Finer grid and smaller material particles will lead to more accurate results. In some simulations, high stress gradients exist in small regions. For instance, in indentation with a sharp tip, the stress gradient is high in the workpiece beneath the indenter tip. In simulation of fracture problems, the stress gradient at the crack tip is high and of particular interest. Consequently, finer grid is needed for these regions; but away from these regions, a coarser grid can be used to reduce the computational cost. In conclusion, a uniform grid can be either too computationally expensive if it is too fine, or inaccurate, if it is coarse. A nonuniform grid with refinement can provide accurate results while minimizing the overall computational time. Grid refinement should maintain the same characteristics of the structured grid as much as possible, in order to replicate the grid generation in parallel processing. The proposed Lx Ly (x1 ,y1) (x2 ,y2) (xp ,yp) 2lx 2ly (x4 ,y4) (x3 ,y3) 45 refinement scheme is illustrated in Fig. 32 with one particle per cell assigned. The material particles that fill each cell are square in nature but denoted as circles for clarity. To understand this grid, one can consider that there are two overlapped structured grids. The coarse grid covers a rectangular region from (0, 0) to (8, 6) and the fine grid covers a region from (2, 2) to (6, 6). For each grid, the shape function can be evaluated from Eq. (312). To be consistent with any other general refinement, it is required that the region of the fine grid to be smaller than the coarse grid. Fig. 32: Refinement of structured grid with a refinement ratio of two When these two grids are merged into one, the shape function and the weighting function for the nodes at the boundary of the finer grid, for example, the nodes at (2, 2) and (2, 3) should be changed. These nodes are called transition nodes. To facilitate the computation of the interpolation function, define an influence zone for each node, denoted as [ − , + , − , + ] x x y y L L L L in 2D or [ − , + , − , + , − , + ] x x y y z z L L L L L L in 3D. The symbols in the square bracket define the size of the influence zone, whereas the subscript denotes the coordinate axis and the superscript denotes the direction of the axis. For example, − x L and + x L represent the sizes in the negative and positive X direction, respectively. The influence zone for each node in 2D is rectangular and it extends to the next immediate grid line to the left, X Y 0 2 4 6 8 4 2 6 46 right, bottom and top of the node. If no more grid lines exist in any direction, such as the boundary nodes, the size is zero in that direction. For example, in the refined grid in Fig. 32, the influence zones for the nodes at (2, 3) and (2, 4) are both [2, 1, 1, 1]. The influence zone for the node at (2, 2) is [2, 1, 2, 1]. Based on this definition, the influence zone for this node in the coarse grid is [2, 2, 2, 2], and in the uniformly fine grid is [0, 1, 0, 1]. Based on the influence zone, the nodal tent function in each direction can be modified as, for example, in the X direction, ⎪ ⎪ ⎩ ⎪ ⎪ ⎨ ⎧ ≤ − − − ≤ − ≤ + − − ≤ − ≤ − ≤ − = + + + − − − x i i x i x i x x i i x x i L x x x x L x x L x x L L x x x x L S x 0 1 ( ) / 0 1 ( ) / 0 0 ( ) . (314) Eq. (314) can be substituted into the grid shape function (Eq. (311)), and the weighting function between the particle p and the node i can be evaluated as ⎪ ⎪ ⎪ ⎪ ⎩ ⎪ ⎪ ⎪ ⎪ ⎨ ⎧ − − − ≥ − − − ≤ − + − ≤ − ≥ = − + + − − + otherwise l b a a L b L a l b a b a L b l b a b a L B L or A L S p x x p x p x x x x ip 2 /(2 ) /(2 ) 0 2 ( ) /(2 ) 0 2 ( ) /(2 ) 0 2 2 2 2 2 2 . (315) where i p A = x − x −l , i p B = x − x + l , a =max( , − − ) x A L and b =min( , + ) x B L . When − = +x x L L , Eqs. (314) and (15) are degraded to the cases for uniform grid. Without detailed proof, the grid shape function and the weighting function still satisfy partition of unity. Similarly, the gradient of the modified weighting function can be computed. 47 It may be noted that in Fig. 32 the refinement ratio is two, i.e., the length of a side of a coarse cell is twice that of the fine cell. To maintain the convenience of the structured grid, only integer refinement ratio should be used. All nodal positions can be computed from the domain of each grid and the cell sizes. The proposed refinement scheme can be applied to any integer refinement ratio and for multiple times for successive refinements. As an example, the weighting function between a particle of size 0.5×0.5 and a node at (0, 0) with an influence zone of [1, 1, 1, 1] is shown in Fig. 33 (a). The particle is on the XY plane and the weighting function is computed at each particle position. For comparison, the influence zone is changed to [1, 0.5, 1, 1], representing a transitional node, while other conditions are the same. The weighting function for this case is plotted in Fig. 33 (b). It can be seen that the weighting function for the transitional node is still C1 continuous. Fig. 33: Effect of influence zone on the weighting functions 3.4 Numerical Simulations 3.4.1 Tracking particle deformation Prior to presenting the results of structured mesh refinement, tracking particle deformation is addressed first since this is necessary in later simulations to achieve accuracy. The material particles are initialized into regular shapes, normally square and (a) Influence zone: [1, 1, 1, 1] (b) Influence zone: [1, 0.5, 1, 1] 1.5 1 0.5 0 0.5 1 1.5 1 0.5 0 0.5 1 1.5 0 0.2 0.4 0.6 0.8 Y Weight X 1.5 1 0.5 0 0.5 1 1.5 1 0.5 0 0.5 1 1.5 0 0.2 0.4 0.6 0.8 Y Weight X 48 cube for 2D and 3D simulations, respectively. All the physical quantities in a particle domain are considered to be uniform. The shape of the particle changes during deformation. So, it is important to track the deformation of each particle. Fig. 34 illustrates the deformation of the particles in 2D when the particles are stretched in the Xdirection. If the particle deformation is not tracked, gaps will form between neighboring particles in the Xdirection, as shown in Fig. 34 (a). Due to the Poisson’s effect, there will be overlapping between particles in the Ydirection, if the particles do not follow the deformation of the materials properly. When the stretch and gaps are large enough, the particles would be separated. Fig. 34 (b) shows the correct deformation in which contiguous particles remain contiguous after deformation. Fig. 34: Schematics showing overlaps and gaps that may occur when particle deformation is not tracked To track the particle deformation, a convenient approach is to calculate the deformed particle shape based on strain history. Since the linear strain increment is computed in GIMP, the effectiveness and validity of tracking particle deformation based on strain would be limited to relatively small deformations. As an example, Fig. 35 (a) shows the simulation of a uniaxial tension problem under plane strain conditions. The material is silicon. Its mass density is 2.71 g/cm3, Young’s modulus 175.8 GPa, and the Poisson’s Overlap X Y (a) No particle deformation (b) Actual particle deformation Gap 49 ratio 0.28. The background grid size is 0.002 × 0.002 mm2. One particle per cell is assigned initially and the time step is 0.02 ns (1 ns = 109 s). The applied pressure increases linearly with time from 0 to 10 ns and is then maintained constant, as shown in Fig. 35 (b). The elongation in the Xdirection of the particle is computed as (1 ) 0x x +ε l , where 0 x l is the initial length of the particles. Separation of the particles occurred at ~25% strain at 6 ns before the full pressure was applied, as shown in Fig. 35 (c). Similar problems have been reported by Guilkey (2005). Fig. 35: Simulation showing separation when the particle deformation is tracked by strain Tracking particle deformation by strain could be more effective if nonlinear strain is used in the GIMP method. However, recovery of the deformed shape based on nonlinear strain involves additional complications. Based on the GIMP algorithm, there is another convenient approach to track the particle deformation. Numerically, the displacement and velocity of each particle in GIMP are computed at the center of the particle to represent X Y L=0.06 mm H=0.04 (a) (c) p (b) t p 40 GPa 0 10 ns X Y 0 0.02 0.04 0.06 0.08 0.1 0.01 0.02 0.03 0.04 0.05 σx 95000 85500 76000 66500 57000 47500 38000 28500 19000 9500 0 Separation (MPa) 50 the entire particle domain. However, in reality, the velocity and deformation at the corners of a particle can be different from the center. Fig. 36 shows four 2D contiguous particles sharing one common corner point at the middle. This corner point should have unique displacement and velocity. As a result, it is helpful to track the displacement and velocity of each corner to track the particle deformation. It is not difficult to compute the velocity of the particle corner given its location. It is computed from the interpolation from the background grid, similar to the center of the particle. For a 2D particle, in addition to updating the position of the center of the particle, the positions of the four corners are updated at each increment. To compute the weighting function for the corners, a fictitious size can be assigned to each corner. Numerical simulation shows that the result is not sensitive to this size in the range of 10% to 80% of the initial particle size. The new particle shape can be obtained by connecting the four corners with straight lines. In order to avoid numerical integration of the interpolation function, it is assumed that the deformed material particle shape is rectangular with edges parallel to the coordinate axes. The size of the rectangle is, therefore, determined from the extent of the corners. As will be demonstrated later in this section, this assumption does not introduce any significant error while it can greatly improve the efficiency of the GIMP algorithm. Fig. 36: Velocity field of a continuum region (the arrows represent both direction and magnitude) 1 2 4 3 51 Fig. 37: GIMP results with tracking deformation of corners and their comparison with FEM Using this approach to track the particle deformation, the problem in Fig. 35 (a) was simulated again with GIMP and the results at 20 ns are plotted in Fig. 37 (a). No separation of particles was seen during the entire simulation up to 50% strain. It is noted that each material particle is plotted as a square of the same size and the particle deformation is not shown due to software limitations on visualization. For comparison, the same problem is simulated using FEM (Abaqus/Explicit) and the FE result is shown in Fig. 37 (b). It can be seen that these two sets of results agree reasonably well with each other; the maximum difference in maximum tensile stress is ~8%. 3.4.2 Indentation problem To verify the refinement algorithm, a 2D indentation problem was simulated. Pressure is applied at the top of the workpiece, as shown in Fig. 38. The dashed lines indicate the (b) FEM (MPa) (a) GIMP 52 borders of refinement levels. The work material is silicon with the properties the same as those given in the previous section. Due to symmetry with respect to the Yaxis, only half of the workpiece is modeled. The size of the model is 0.027 × 0.04 mm2. The pressure increases linearly with time from t = 0  30 ns and is then kept constant at 60 GPa. Fig. 38: Twodimensional indentation simulation Several simulations were performed under different settings for the purpose of comparison. In the first simulation, a uniform grid with a cell size of 0.001 × 0.001 mm2 is used and the time step is 20 ps (1 ps = 10–12 s). The stress distribution in the workpiece at t = 20 ns is shown in Fig. 39 (a). In this figure, the units of length and stress are mm and MPa, respectively. In the second simulation, as indicated in Fig. 39 (b), two levels of refinements are used with the refinement ratio to be 2. The cell lengths are 0.002 mm and 0.001 mm for the first and second levels, respectively. The fine level covers a rectangular area of the workpiece from (0, 0.02) to (0.02, 0.04). The grid is fixed in space; the material particles initially in the fine region move to the coarse region during deformation. As shown in Fig. 39 (b), some fine material particles have moved below the line Y = 0.02 mm. It may be noted that in Fig. 39 the material particles are plotted as squares corresponding to their initial sizes. Gaps between particles are intentionally shown to depict the particle sizes. X Y p Level 2 Level 3 Symmetry p t 60 GPa 0 30 ns Level 1 53 Fig. 39: Comparison of the stress distributions at different levels of refinements in force indentation Three levels of refinements are used in the third simulation of the same problem. In this simulation, the time step is 10 ps and the results are shown in Fig. 39 (c). The stresses agree very well with the previous two simulations. Fig. 39 (d) shows the result when the particle deformation is not tracked. A severe material separation was observed at t = 36 ns, as indicated by the arrow. Additionally, the displacement history of the particle in the middle of the top surface for the third simulations is shown in Fig. 310. The same result of the integration point of the element in the middle of the top surface in FE simulation using ABAQUS/Explicit is also shown in Fig. 310 for comparison. It can be seen that Separation (a) One level (b) Two levels (c) Three levels (d) Deformation not tracked (MPa) 54 the displacement compares well with FE up to t = ~20 ns. After this time, FE simulation aborted due to mesh distortion. This demonstrates the capability of GIMP using a grid with structured refinement in handling large deformations. Time (ns) Uy (μm) 0 10 20 30 40 8 6 4 2 0 GIMP FEM FEM aborted Fig. 310: Comparison of displacement history with FE It may be noted that if the exterior corners of the surface particles are tracked from the nodal interpolations in the same way as the interior corners, simulations tend to become unstable due to erroneous surface corner displacements. This problem was observed to be strictly and consistently associated with the particles with external tractions applied. It is caused by insufficient nodal interpolation. To eliminate this problem, the exterior corners of the surface particles were tracked by strain only, as used in these simulations. 3.4.3 Verification of the displacement boundary condition Next, the results on the validation of the displacement boundary conditions are presented. For dynamic simulations, an artificial damping may be introduced. With damping, the nodal momentum can be updated as d t i u i i b i i i Δp = (f int + f + f τ + f − p )Δ , (316) where d is the artificial damping coefficient. 55 A rectangular slab is fixed on the left and a displacement boundary condition is applied on the right, as shown in Fig. 311. The material is silicon and its properties are given in the previous section. The cell size in this simulation is 0.5 × 0.5 mm2 and four particles are assigned to each cell initially. The time step is 5 ns. The prescribed displacement increases linearly with time to 0.5 mm at t = 10 μs, which corresponds to a velocity of 50 m/s, and remains constant thereafter. This problem is simulated in FE using Abaqus/Explicit for comparison. The displacement in the Xdirection, Ux, for a particle initially centered at (2.625, 1.625) as a function of time is plotted in Fig. 312 (a). It can be seen that without damping, the vibrations of displacement from both FEM and GIMP simulations are in phase before t = 14 μs, but out of phase afterwards. The dashed line represents the steady state displacement at this point. It can be seen that when an artificial damping of 106 s1 is used, the GIMP solution converges quickly. The error of the converged displacement is 0.46%. Fig. 312 (b) shows the comparison of the normal stress in the Xdirection. Good comparison between GIMP and FE has been obtained. With damping, the GIMP solution converges to the analytical solution for static simulations. Fig. 311: A simple tension problem with displacement boundary conditions X Y L=6 mm H=3 mm u 56 Fig. 312: Comparisons of the displacement and stress 3.4.4 Stress concentration problem Fig. 313 shows a copper plate (60 × 60 mm2) with a central hole (4 mm diameter) subjected to a distributed load. This problem is simulated using GIMP as a dynamic problem with a damping factor of 1000 s1 and three levels of refinement. The cell sizes at these three levels are 1.0 mm, 0.5 mm, and 0.25 mm, respectively. One particle is assigned to each cell not adjacent to the hole initially. The cells close to the circular hole are assigned 25 particles each to model the circular edge more accurately with the use of square particles. It may be noted that all the particles occupy square areas initially. The time step is 10 ns and the applied distributed load intensity p = 10 MPa. Fig. 313: A plate with a circular hole subjected to tension The stress distribution after 4000 steps is shown in Fig. 314 (a) and the area close to the hole (bottom left corner in Fig. 314 (a)) is magnified in Fig. 314 (b). The normal stress p p X Y θ Time (μs) Ux (mm) 0 10 20 30 40 0 0.05 0.1 0.15 0.2 0.25 0.3 GIMP d=0 FE d=0 GIMP d=106 1/s Steady State (a) Displacement history Time (μs) σx 0 10 20 30 40 0 5 10 15 20 GIMP d=0 FE d=0 GIMP d=106 1/s Steady State (GPa) (b) Stress history 57 of the particle at the top of the hole is 30.7 MPa when the applied tension is 10 MPa. This gives a stress concentration factor of 3.07. The tangential stress of the particles along the hole circumference, normalized by the applied pressure, is plotted in Fig. 315 in comparison with the theoretical value (Pilkey (1997)). A good agreement between the numerical simulation and theoretical value is obtained. This demonstrates that the GIMP refinement algorithm is effective for problems involving significant stress gradients. Furthermore, with the use of small square particles, GIMP is capable of modeling curved surfaces . Fig. 314: Normal stress in the Xdirection with p = 10 MPa Angle θ (degrees) Normalized Tangential Stress 20 15 30 45 60 75 90 1 0 1 2 3 4 Numerical Theoretical (Pilkey, 1997) Fig. 315: Normalized tangential stress along the circumference of the hole (a) Overall (b) Magnified (MPa) 58 3.4.5 Static stress intensity factor Next, the GIMP refinement algorithm is implemented to determine the stress field and stress intensity factor for a mode I crack problem to determine the capability of the GIMP refinement algorithm in simulation of stress distribution near a crack. Guo and Nairn (2003) have successfully extended the MPM method to compute stress distribution in a plate with explicit cracks using multiple nodal fields along the crack surface. The physical quantities of material particles on each side of the crack were interpolated using variables in the field on that side of crack surface. In their simulation, a uniform mesh was used. Since the stress gradient at the crack tip is very high, a refined mesh near the crack tip and a coarse mesh in the far field should lead to savings in computational time while maintaining the same accuracy with the use of a uniform fine mesh. The fracture problem is thus an appropriate problem to evaluate the refined GIMP algorithm. For this purpose, the same fracture problem by Guo and Nairn (2003) using MPM is modeled, that is, a double cantilever beam (DCB) with a crack as shown in Fig. 316. Fig. 316: Geometry of a double cantilever beam with a crack In the area close to the crack tip, finer meshes are used, while in the area far away from the crack tip, coarse meshes are used. The thickness of the plate is 1 mm, thereby justifying a plane stress condition for this problem. The material of the DCB is L=10 2H=2 a=5 X Y F F Unit: mm Crack 59 considered to be homogeneous, isotropic, and linearly elastic. Its density, Young’s modulus, and Poisson’s ratio are 1500 kg/m3, 2300 GPa, and 0.33, respectively. The applied force is F = 4×10−4 N and this results in a mode I crack problem. The static stress intensity factor for the DCB can be calculated using the following equation (Kanninen 1973) 3 / 2 2 3 ( 3 / 2) H K F a H I + = . (317) This problem has been simulated using MPM with uniform grids of three sizes, 4 mm, 2 mm, and 1 mm (Guo and Nairn (2003)). They have computed the energy release rate and determined the stress intensity factor from Jintegral. Their results indicate that the stress intensity factor determined from finer grid is more accurate and closer to the theoretical value. In the mesh refinement GIMP algorithm used in this study, the energy release rate was computed using the virtual crack closure technique (Rybicki and Kanninen (1977); Wang, Karuppiah, Lu, Roy, and Komanduri (2005)). The energy released during an infinitesimal crack growth of Δa is assumed to be the energy required to close the crack to its initial size. Hence, the energy release rate G is determined by ∫Δ Δ → ⋅ Δ Δ = a a da a G 0 0 ( ) 2 lim 1 σn u , (318) where n is the crack surface normal. For the 2D mode I crack in the Xdirection, ∫Δ Δ → ⋅Δ Δ = a I a yy yσ u da a G 0 2 0 lim 1 . (319) In GIMP, the energy release rate can be computed as 60 ( ) 2 1 1 2 I tip y y F u u t a G − Δ = , (320) Fig. 317: Material points and background grid around the crack tip where the superscripts 1 and 2 denote the two material particles immediately to the left of the crack tip as shown in Fig. 317, Δa is the X distance between particle 1 and the crack tip, t is the thickness of the beam. tip F is the nodal force to hold the crack tip together (Rybicki and Kanninen (1977)) and is computed as the crack tip nodal force from one side of the crack in this simulation. The mode I stress intensity factor for the static crack is given by K GE I = . (321) GIMP simulations were carried out using four material points per cell. The time step is 0.1 μs, and a damping coefficient of 4000 s1 is used to allow the results to converge to the static data. The computed stress intensity factors using uniform grids of 4 mm and 1 mm, respectively, are plotted in Fig. 318 (a) and the theoretical value calculated from Eq. (317) is also shown for comparison. For simulation with refinement, two levels of refinement were used, i.e., 2 mm grid size for the coarse level and 1 mm grid size for the fine level, and the extent of the fine level is indicated by the dashed square in Fig. 316. The computed stress intensity factor using the refinement algorithm is also plotted in Fig. 318 (a). It is seen that the results from two levels of refinement are identical to the Crack surface X Y 1 2 61 results from one uniform fine level with 1 mm cell size. Moreover, the simulation time using the refinement algorithm is 38% shorter than that of one uniform fine level. The computed stress intensity factor became even closer to the theoretical value when a third level of refinement was added for the crack tip as shown in Fig. 318 (b). For the case of using uniform grid of 2 mm, the computed stress intensity factor using GIMP and MPM (Guo and Nairn (2003)) are plotted in Fig. 319. The MPM results were computed using a damping factor of 1000 s1, and therefore, more oscillations can be seen as expected. Fig. 318: Computed stress intensity factor as a function of time Fig. 319: Comparison of the stress intensity factor with MPM results (a) Overall (b) Magnified 62 Fig. 320: Stress distribution in the beam using three levels of refinements at t = 2.5 ms Fig. 320 shows the distribution of σy at t = 2.5 ms when the force applied on the beam was changed to F = 4 N, 10000 times of the previous value, with other parameters the same. In this simulation, three levels of refinements with cell sizes of 1 mm, 0.5 mm and 0.25 mm, were used and the different density of material points in each level can be clearly seen in the figure. The computed stress intensity factor, scaled by 10000 times, still compares well with the theoretical value. It is noted that deformation near the crack surfaces has caused material points crossing cell boundaries, a situation where MPM would give numerical noise such as alternating stress signs. Despite the extent of the deformation, and material points crossing cell boundaries, GIMP method with the use of particle deformation tracking still gives correct results, further verifying the structured refinement methods developed herein. 3.5 Conclusions 1. A spatial refinement scheme for a structured grid was developed by adding transitional nodes and by changing the influence zone of the transition nodes in GIMP. The influence zone is square for uniform grid nodes and rectangular for transitional nodes. This influence zone affects the computation of the nodal shape function. The computation of the weighting function remains the same as for the uniform grid. The refinement scheme (MPa) Crack 63 can be applied successively and the refined grid remains structured in each refinement level, i.e., every node can be determined by the extent of the grid level and cell size. The refinement scheme was implemented and several problems such as tension, indentation, stress concentration, and stress distribution near a crack (mode I crack problem) were modeled to demonstrate its effectiveness and accuracy. A good agreement has been obtained between numerical and theoretical results, indicating the validity of the structured mesh refinement for GIMP scheme. 2. The GIMP algorithm has also been extended to include the displacement boundary condition, based on the approach used in the meshless local PetrovGalerkin (MLPG) method, Atluri and Zhu (2000). A penalty parameter is used to impose the displacement boundary condition and a nodal force vector because the displacement boundary condition is introduced to the nodal momentum governing equation. A uniaxial tension problem with constant pulling velocity was simulated to verify the displacement boundary condition. 3. A method to track the material particle deformation was developed and verified in one example. When the particle deformation is not tracked, artificial separation was observed when the particle strain increases to a certain level. In tensile simulations, when the normal strain is ~25%, material particles tend to separate from the body. Our approach tracks the displacement of each corner of the material particles. Since neighboring particles share corners, no separation would occur during deformation using this approach. 64 Chapter 4 Simulation of Deformation Mechanism and Failure of Bulk Metallic Glass (BMG) Foam Using the GIMP Method A new bulk metallic glass (BMG) foam, processed by thermoplastic expansion, has recently been discovered (Schroers et al. (2003); Brothers and Dunand (2004); Hanan et al. (2005)). In this investigation, the microstructure and mechanical properties of BMG foam were determined using microtomography combined with insitu compression. Numerical simulation of the mechanisms of deformation and failure of the BMG foam, which has never been done before, is conducted in this investigation. The complex cellular geometry of the foam was modeled using the generalized interpolation material point (GIMP) method by creating material points based on grayscales at voxels determined from microtomography. Nanoindentation tests were conducted to determine the local elastic modulus of the foam matrix and the results are used as the input to the GIMP simulation. Several cubic volume elements (from 0.5763 mm3 to 1.7283 mm3) taken from the 3D tomographic image were used in the simulation. The size of the cube is increased until further increase in volume does not lead to a change in mechanical response; the cube with this size is determined as the representative volume element. The simulated compressive stressstrain curve using GIMP was compared with the experimental data and a good agreement on the initial slope (representative of the average Young’s modulus of the foam) was obtained. The stress contour in the foam was 65 illustrated by analyzing different sections of the foam at different simulation times. The densification and hardening of the ductile metallic foam under compression was simulated to demonstrate the capabilities of GIMP and to study the failure mechanisms. 4.1 Introduction Metallic glass foam has been used in various engineering applications due to such combination of unique properties as low density, high specific strength, high energy absorption, and superior thermal insulation. Foam materials are generally divided into opencell and closedcell foams. In general, opencell foams have higher porosity and softer behavior than closedcell foams. The uniaxial mechanical behavior of a foam is characterized by three regimes, namely, initial elastic response, the stress plateau that follows, and final stiffening due to densification (Gibson and Ashby (1997); Miller (2000); Katti et al. (2006)). Theoretical modeling of the foam mechanics has been carried out for both regular and irregular opencell foams. The nonlinear stressstrain relation of the foams at high strains was modeled by the buckling of elastic cell edges (Zhu et al. (1997)) or by the implementation of nonlinear kinematics (Wang and Cuitino (2000)). Most numerical simulations of foam materials reported in the literature were either at the macroscopic scale or the mesomechanical (micromechanical) scale. At the macroscopic scale, the foam material is considered as a homogeneous material with homogenized properties determined from measurements such as uniaxial compressions (Gibson (1997)). The work using this approach is rather limited, due primarily to the wellknown problem associated with simulating cells in contact. In this approach, large foam structures can be simulated at relatively low computational cost and the deformations under different loading conditions can be predicted reasonably well (Meguid et al. (2002); Meguid et al. 66 (2004)). However, such simulations require knowledge of loading conditions a priori and the associated material model may not be readily available in some situations (Wicklein et al. (2004)). In the mesomechanical simulations, the microstructure of the foam is modeled explicitly. Xray tomography is used to determine the threedimensional (3D) microstructure and the 3D geometry is used in a numerical model for simulation. In this approach, the discrete element/particle size is dictated by the resolution of the Xray tomography. Kadar et al. (2004) investigated the mechanical behavior of closedcell Al porous foam by indentation using the finite element method (FEM). Wicklein and Thoma (2005) conducted a series of FEM simulations and determined the elastic and plastic response of an opencell aluminum foam. The mesomechanical approach provides microstructural evolution during deformation; however, the computational cost can be very high. For model generation using mesomechanical simulation approach, Wicklein et al. (2004) described three possible approaches to discretize a foam structure into FEM meshes. In the first approach, beam and plate elements are constructed based on the cell votices in the tomographic images. In the second approach, voxels are transformed into cubic elements in a predefined structured mesh. In the third approach, tetrahedral elements are used to mesh the entire foam structure. The second approach was used successfully for the simulation of aluminum foams (Wicklein et al. (2004); Wicklein and Thoma (2005)) under a few percent deformations. However, simulation of a fully densified foam is a challenge in FEM because of large distortions of the microstructure involved and the internal contact of cell walls upon closure of cells. Brown et al. (2000) used hexahedral elements to model the struts of an opencell foam under impact loading. It was pointed 67 out that the drawbacks of this approach involve the use of a large number of elements and small timesteps. It is noted that the struts all have regular cross sections so that it was relatively easy to mesh all the struts. Recently, the generalized interpolation material point (GIMP) method was introduced for dynamic simulation of materials (Bardenhagen and Kober (2004); Ma et al. (2005)). GIMP has evolved from the material point method (MPM), originally developed by Sulsky, Zhou and Schreyer (1995) and subsequently refined by others (e.g., Tan and Nairn (2002); Wang et al. (2005)). In GIMP, the workpiece is discretized into a collection of material particles (points). The material points carry all the physical variables, including the mass, momentum and stress, in the simulation. A background grid, usually a structured grid fixed in space, is used to discretize the momentum conservation equation and solve it at the grid nodes using explicit time integration. In a regular GIMP computational step, the material point information is first interpolated to the grid nodes. Then, the momentum conservation equation is solved at the nodes to update the nodal information. Next, the updated nodal information is interpolated to update the position, velocity, strain, and stress at the material points. Finally, the background grid is redefined (usually reset the nodal quantities to zero but keep the nodal positions fixed). GIMP shares some similarities with other meshless methods, e.g., Atluri and Zhu (2000). Fig. 41 (a) illustrates an elliptical workpiece discretized into a collection of material points of different sizes with the background grid shown. When two bodies are close enough, as indicated by open and closed circles in Fig. 41 (b), the material point information will be interpolated to the same grid node in the middle when one background is used. In this case, these two workpieces are in nonslip contact. The 68 intrinsic features associated with GIMP, such as ease of discretization, no bodyfixed meshes used so that mesh distortion is not an issue and natural capability of handling nonslip contact, made the simulation of some complex problems possible. Fig. 41: Illustrations of the discretization scheme and intrinsic contact in GIMP GIMP was used successfully to simulate densification of opencelled foam materials (Bardenhagen et al. (2005); Brydon et al. (2005)). In GIMP simulations, each voxel in the Xray tomograph is converted into a material particle and an Eulerian structured grid is used for solving the field equations. A representative volume element (RVE) was used to simulate the bulk response of a foam. Their results have indicated that the apparent Poisson’s ratio is negative in the stress plateau and the dynamic response of the compression is dominated by a compression wave, which is an inertial effect absent in homogeneous material and its velocity is much slower than any characteristic wave speeds (Brydon et al. (2005)). In this work, we report the results of a simulation of a closedcell foam made of bulk metallic glass. Most metal and its alloys in engineering applications are in the form of crystalline structures; their mechanical properties are dependent on the crystal structure, chemistry, microstructure, and processing conditions. Metallic glass is another form of a metal or an alloy. In a metallic glass, the material is amorphous and has no longrange atomic order. Fig. 42 (a) and (b) are TEM images of the microstructure of a crystalline (a) Discretization and background grid (b) Nonslip contact 69 zirconium alloy and the same alloy in the amorphous form (Hufnagel (2005)). Metallic glass can be formed by rapid cooling of a liquid metal; the cooling rate is so high that ordered crystalline structure cannot be developed in the metal, resulting in amorphous structures. Fig. 42: High resolution TEM images of (a) crystalline and (b) amorphous structures of a zirconium alloy (Hufnagel (2005)) Metallic glass foams recently have attracted great attention due to their unique properties. They exhibit high specific strength, high stiffness, high energy absorption (Ashby et al. (2000)), and more importantly they can be used under those conditions, including, high temperatures where other foams (e.g., polymer foams) cannot be used. BMG foams combine exceptional strength, elasticity, wear, and corrosion resistance with modest densities and low processing temperatures. They are being considered for structural and biocompatible implant applications (Hanan et al. (2005); Brother and Dunand (2005)). The thickness of BMG material is usually limited (on the order of a few millimeters to a centimeter or so) by the requirement that a high enough cooling rate in the center of the wall is necessary to form the amorphous structures. For a
Click tabs to swap between content that is broken into logical sections.
Rating  
Title  Multiscale Simulation Using the Generalized Interpolation Material Point Method, Discrete Dislocations and Molecular Dynamics 
Date  20060501 
Author  Ma, Jin 
Department  Mechanical Engineering 
Document Type  
Full Text Type  Open Access 
Abstract  A multiscale simulation scheme that spans the atomistic scale to the continuum has been developed for materials simulations in this study. At the continuum scale, the generalized interpolation material point (GIMP) method has been extended for parallel processing using the Structured Adaptive Mesh Refinement Application Infrastructure (SAMRAI). A contact algorithm in GIMP has been developed for the treatment of contact pair between a rigid indenter and a deformable workpiece. Two spatial refinement schemes for GIMP are presented for simulations with highly localized stress gradients at the continuum scale. A method for multiscale simulation bridging different scales, namely the continuum scale using GIMP, the mesoscale using discrete dislocations and the atomistic scale using the molecular dynamics (MD), is presented and validated in two dimensions. Findings and conclusion/. Numerical simulations with multiple length scales from nanometer to millimeter were conducted and validated on a 2D nanoindentation problem. Numerical simulations of several problems, such as tension, indentation, stress concentration and stress distribution near a crack (mode I crack problem) are presented to validate the refinement schemes at the continuum scale as well as the parallel processing algorithm. The capability of handling large deformation in GIMP is also demonstrated. A mode I crack propagation problem is simulated using the coupling algorithm. The stress field near the crack tip was validated by comparing results from coupled simulations with purely GIMP simulations of the same model. Coupled simulation results were also compared with purely MD simulations. A very good agreement was obtained. Other problems, such as dynamic friction problem at atomistic scale and the nanoindentation problem, are also simulated using the multiscale simulation algorithm. 
Note  Dissertation 
Rights  © Oklahoma Agricultural and Mechanical Board of Regents 
Transcript  MULTISCALE SIMULATION USING THE GENERALIZED INTERPOLATION MATERIAL POINT METHOD, DISCRETE DISLOCATIONS AND MOLECULAR DYNAMICS By JIN MA BACHELOR OF SCIENCE Wuhan University of Technology Wuhan, China June, 2000 MASTER OF SCIENCE Oklahoma State University Stillwater, Oklahoma August, 2002 Submitted to the Faculty of the Graduate College of Oklahoma State University in partial fulfillment of the requirements for the Degree of DOCTOR OF PHILOSOPHY May, 2006 ii MULTISCALE SIMULATION USING THE GENERALIZED INTERPOLATION MATERIAL POINT METHOD, DISCRETE DISLOCATIONS AND MOLECULAR DYNAMICS Dean of the Graduate College Thesis Advisor Dissertation Approved: Hongbing Lu, Ranga Komanduri Jay C. Hanan Lionel M. Raff Demirkan Coker A. Gordon Emslie iii © Copyright By JIN MA May 05, 2006 iv Acknowledgements I wish to express my most sincere appreciation to my thesis advisor, Dr. Hongbing Lu, for his intelligent supervision, constructive guidance, inspiration and friendship. Dr. Lu has also given me very generous financial support. My sincere appreciation extends to other committee members, Dr. Ranga Komanduri, Dr. Demir Coker and Dr. Lionel Raff, whose guidance, assistance, encouragement, and friendship are also invaluable. I would like to thank Dr. Bo Wang, Dr. Jay Hanan, Dr. Samit Roy, Dr. Martin Hagan, Dr. Afshin J. Ghajar, Dr. C. Eric Price and Dr. J. Keith Good for their long time assistance and guidance. Moreover, I wish to express my sincere gratitude to those provided suggestions and assistance for this study: Mr. Gang Huang, Ms. Haowen Yu, Mr. Nitin Daphalapurkar and Dr. Huiyang Luo. My research work was supported by a grant from the Air Force Office of Scientific Research (AFOSR) through a DEPSCoR grant (No. F496200310281). I would like to thank Dr. Craig S. Hartley, and Dr. Jamie Tiley, Program Managers for the Metallic Materials Program at AFOSR for the support of and interest in this work. Thanks are also due to Dr. Richard Hornung and Dr. Andy Wissink of the Lawrence Livermore National Laboratory for providing the code and assistance with the use of SAMRAI, Dr. Scott Bardenhagen for his introduction to the GIMP algorithm and valuable comments, and Dr. John A. Nairn for providing the 2D MPM code. I want to give my special appreciation to my wife, Yang Liu, for her precious suggestions v to my study, her encouragement at times of difficulty, her love and understanding throughout my study. Thanks also go to my parents and my friends for their support and encouragement. Finally, I would like to thank the School of Mechanical and Aerospace Engineering for supporting me during my M.S. and Ph.D. studies. vi Table of Contents Acknowledgements............................................................................................................ iv Chapter 1 Introduction ........................................................................................................ 1 Chapter 2 Multiscale Simulations Using Generalized Interpolation Material Point (GIMP) Method and SAMRAI Parallel Processing ......................................................................... 4 2.1 Introduction.........................................................................................................5 2.2 Generalized Interpolation Material Point (GIMP) Method ................................ 9 2.2.1 Contact algorithm in GIMP ...................................................................... 12 2.2.2 Numerical implementation........................................................................16 2.3 Parallel Computing Scheme Using GIMP with SAMRAI ............................... 17 2.3.1 SAMRAI...................................................................................................17 2.3.2 Spatial and temporal refinements.............................................................. 18 2.3.3 Domain decomposition .............................................................................22 2.4 Numerical Examples and Results ..................................................................... 25 2.5 Conclusions.......................................................................................................34 Chapter 3 Structured Mesh Refinement in GIMP Method for Simulation of Dynamic Problems ........................................................................................................................... 37 3.1 Introduction.......................................................................................................38 3.2 GIMP.................................................................................................................40 3.3 Structured Mesh Refinement ............................................................................44 3.4 Numerical Simulations......................................................................................47 3.4.1 Tracking particle deformation................................................................... 47 3.4.2 Indentation problem..................................................................................51 3.4.3 Verification of the displacement boundary condition............................... 54 3.4.4 Stress concentration problem.................................................................... 56 3.4.5 Static stress intensity factor ...................................................................... 58 3.5 Conclusions.......................................................................................................62 Chapter 4 Simulation of Deformation Mechanism and Failure of Bulk Metallic Glass (BMG) Foam Using the GIMP Method............................................................................ 64 4.1 Introduction.......................................................................................................65 4.2 Nanoindentation Testing of BMG Foam .......................................................... 70 4.3 Reconstruction of the BMG Foam and Simulation Convergence .................... 72 vii 4.3.1 Reconstruction of the BMG foam for GIMP ............................................ 72 4.3.2 Effect of grid cell size on simulation convergence................................... 74 4.4 Simulation of the BMG Foam in Compression ................................................ 76 4.4.1 Material failure..........................................................................................78 4.4.2 Stress flow in the foam ............................................................................. 79 4.4.3 The minimum size of a representative volume element (RVE)................ 81 4.4.4 Simulation of densification ....................................................................... 84 4.5 Conclusions.......................................................................................................89 Chapter 5 Multiscale Simulation Using GIMP Method and Molecular Dynamics .......... 91 5.1 Introduction.......................................................................................................92 5.2 GIMP and Refinement ...................................................................................... 96 5.3 MD Simulation and Atomistic Strain ............................................................... 99 5.3.1 Incremental atomic strain........................................................................100 5.3.2 Interpolation function..............................................................................101 5.3.3 Numerical validation...............................................................................104 5.4 Coupling of GIMP and MD Simulations ........................................................ 107 5.4.1 Coupling scheme.....................................................................................107 5.4.2 Parallel processing of the coupled model ............................................... 111 5.5 Numerical Results........................................................................................... 112 5.5.1 Multiscale simulation of mode I crack.................................................... 112 5.5.2 Multiscale simulation of mode II crack .................................................. 117 5.6 Conclusions..................................................................................................... 119 Chapter 6 Multiscale Simulation of Nanoindentation using GIMP, Dislocations Dynamics and Molecular Dynamics ................................................................................................ 122 6.1 Introduction..................................................................................................... 122 6.2 Coupling Scheme between GIMP, DD and MD............................................. 126 6.2.1 Coupling of GIMP and MD .................................................................... 126 6.2.2 Bridging the continuum and atomistic scales with DD .......................... 128 6.2.3 Parallel processing ..................................................................................130 6.2.4 Contact solver for MD ............................................................................ 133 6.3 Simulation of Nanoindentation ....................................................................... 133 6.4 Conclusions..................................................................................................... 140 Chapter 7 Summary and Future Work............................................................................ 141 7.1 Summary ......................................................................................................... 141 7.2 Future Work.................................................................................................... 143 7.2.1 Other multiscale problems ...................................................................... 143 7.2.2 Simulation of metal cutting using GIMP................................................ 144 viii Disclaimer on SAMRAI ................................................................................................. 149 References...................................................................................................................... 150 ix List of Figures Fig. 21: Illustration of early contact in multimesh ......................................................... 13 Fig. 22: Schematic of contact algorithm between the rigid indenter and the deformable workpiece......................................................................................................................... 14 Fig. 23: Material points in cells and the weighting function in 2D GIMP ...................... 16 Fig. 24: Illustration of a hierarchy of three nested grid levels of mesh refinement......... 18 Fig. 25: Two neighboring coarse and fine grid levels in 2D GIMP computations.......... 21 Fig. 26: Nested multilevel refinement and reduction in the number of cells with number of levels............................................................................................................................ 21 Fig. 27: A computational domain of two patches in one grid level................................. 23 Fig. 28: Flowchart showing advancement of grid levels recursively starting from the coarsest to the finest level in GIMP.................................................................................. 24 Fig. 29: Simulation results of tensile stress contours for a simple tension problem ....... 26 Fig. 210: Loading conditions for a simple indentation problem ..................................... 27 Fig. 211: GIMP and FEM results of normal stress variation in the Ydirection at different increments .......................................................................................................... 28 Fig. 212: Schematic of 2D indentation showing (a) three levels of refinement and (b) the indenter velocity history ................................................................................................... 29 Fig. 213: Comparison of normal stresses in Ydirection from FEM and GIMP............. 30 Fig. 214: Comparison of shear stresses from FEM and GIMP ....................................... 31 Fig. 215: Normal and shear stresses of GIMP simulations with a uniform cell size of 500 nm .............................................................................................................................. 32 Fig. 216: Indentation load versus depth curves from FEM and GIMP with different grid sizes.................................................................................................................................. 33 Fig. 217: Multiscale simulation of nanoindentation with five levels of refinement........ 34 Fig. 31: 2D representation of a particle and a grid cell ................................................... 44 Fig. 32: Refinement of structured grid with a refinement ratio of two ........................... 45 Fig. 33: Effect of influence zone on the weighting functions ......................................... 47 Fig. 34: Schematics showing overlaps and gaps that may occur when particle deformation is not tracked ................................................................................................ 48 Fig. 35: Simulation showing separation when the particle deformation is tracked by strain................................................................................................................................. 49 Fig. 36: Velocity field of a continuum region (the arrows represent both direction and magnitude) ........................................................................................................................ 50 x Fig. 37: GIMP results with tracking deformation of corners and their comparison with FEM ................................................................................................................................. 51 Fig. 38: Twodimensional indentation simulation........................................................... 52 Fig. 39: Comparison of the stress distributions at different levels of refinements in force indentation........................................................................................................................ 53 Fig. 310: Comparison of displacement history with FE.................................................. 54 Fig. 311: A simple tension problem with displacement boundary conditions ................ 55 Fig. 312: Comparisons of the displacement and stress.................................................... 56 Fig. 313: A plate with a circular hole subjected to tension ............................................. 56 Fig. 314: Normal stress in the Xdirection with p = 10 MPa .......................................... 57 Fig. 315: Normalized tangential stress along the circumference of the hole .................. 57 Fig. 316: Geometry of a double cantilever beam with a crack........................................ 58 Fig. 317: Material points and background grid around the crack tip .............................. 60 Fig. 318: Computed stress intensity factor as a function of time .................................... 61 Fig. 319: Comparison of the stress intensity factor with MPM results .......................... 61 Fig. 320: Stress distribution in the beam using three levels of refinements at t = 2.5 ms62 Fig. 41: Illustrations of the discretization scheme and intrinsic contact in GIMP .......... 68 Fig. 42: High resolution TEM images of (a) crystalline and (b) amorphous structures of a zirconium alloy (Hufnagel (2005)) ................................................................................ 69 Fig. 43: Average nanoindentation loaddisplacement curve with error bars for the BMG foam (Pd43Ni10Cu27P20)..................................................................................................... 71 Fig. 44: Schematic of Xray tomography system ............................................................ 72 Fig. 45: 3D microtomographic image of a BMG foam................................................... 73 Fig. 46: Stress distribution on a section simulated with the grid size of 3 and 4 times of the voxel size, respectively ............................................................................................... 75 Fig. 47: Stress distribution on a section simulated with the grid size to be 1 and 2 times of the voxel size. ............................................................................................................... 76 Fig. 48: Experimental setup and nominal stressstrain curve.......................................... 77 Fig. 49: Stress distribution in the overall GIMP model at t = 2.5 μs at location 1.......... 79 Fig. 410: Stress distribution on different sections at t = 2.5 μs. ...................................... 80 Fig. 411: Stress distribution on different sections at 8% strain. ...................................... 80 Fig. 412: Stressstrain curves of the BMG foam from experiments and simulations using second order interpolation................................................................................................. 82 Fig. 413: Stressstrain curves of the BMG foam from experiments, model prediction and simulations using linear interpolation............................................................................... 83 Fig. 414: Stressstrain curve from the densification simulation...................................... 85 Fig. 415: Comparison of initial and final shapes of the RVE.......................................... 86 Fig. 416: Maximum shear stress on section Y50 at different compressive strains ........ 87 Fig. 417: Distribution of z σ stress on section Y50 at different compressive strains .... 88 xi Fig. 418: Histograms of z ε at three compressive strains ................................................ 89 Fig. 51: Two neighboring coarse and fine grid levels in 2D GIMP computations.......... 98 Fig. 52: 2D representation of a particle and a grid cell ................................................. 103 Fig. 53: The generalized interpolation function in 2D .................................................. 104 Fig. 54: Two examples to calculate the atomistic strain: (a) Tension (b) Shear ........... 105 Fig. 55: Strain histories from different simulations....................................................... 106 Fig. 56: Comparison of the normal strain of a surface atom in tension ........................ 106 Fig. 57: Illustration of coupled GIMP and MD simulations. The circles represent atoms and squares (smaller than physical size) material points. The material points connect to each other without a gap to represent continuum. .......................................................... 108 Fig. 58: Flow chart of the coupling algorithm for one increment ................................. 110 Fig. 59: Illustration of three GIMP grid refinements and the domain decomposition for coupling.......................................................................................................................... 112 Fig. 510: Coupled GIMP/MD simulation of a 2D modeI edge crack. The dashed lines are the boundaries of the atomistic domain .................................................................... 112 Fig. 511: Stress distribution for P = 0.3 eV/Å3.............................................................. 114 Fig. 512: Simulation model with an edge crack with P = 0.3 eV/Å3 at t = 208 τ.......... 115 Fig. 513: Energies in the model with P = 0.15 eV/Å3 ................................................... 116 Fig. 514: Energy release rate P = 0.15 eV/Å3................................................................ 116 Fig. 515: Boundary conditions of the mode II crack problem....................................... 117 Fig. 516: Comparison of results from pure MD and coupled simulations at t = 11.2 τ 118 Fig. 517: Dislocation path at different simulation time................................................. 119 Fig. 61: Illustration of coupled GIMP and MD simulations. The circles represent atoms and squares (smaller than physical size) material points. The material points connect to each other without a gap to represent continuum ........................................................... 127 Fig. 62: Illustration of the domain decomposition and refinement for the coupling simulation of 2D indentation using GIMP, DD and MD................................................ 131 Fig. 63: Flow chart of the coupling algorithm for each increment................................ 132 Fig. 64: 2D modeling of an indentation problem with a cylindrical indenter ............... 134 Fig. 65: Dislocations and the loaddepth curve ............................................................. 134 Fig. 66: Comparison of loaddepth curves using different contact algorithms ............. 135 Fig. 67: Effect of the indentation velocity on the loaddepth curve .............................. 136 Fig. 68: Effect of the indenter size on the loaddepth curve.......................................... 137 Fig. 69: Loaddepth curves when the ratio between the workpiece size and the indenter radius is constant............................................................................................................. 138 Fig. 610: Comparison of the loaddepth curves between the coupling and purely MD simulations ...................................................................................................................... 139 Fig. 611: Slip bands, discrete dislocations and stress distributions in the model.......... 139 Fig. 71: Simulation snap shots of 2D orthogonal metal cutting using GIMP ............... 147 xii List of Tables Table 41: Nanoindentation results for the BMG foam. ................................................... 71 Table 42: Porosity of different volume elements at each location. ................................. 81 Table 71: JohnsonCook material parameters for copper (ABAQUS (2005)).............. 146 Table 72: Material properties for copper....................................................................... 146 1 Chapter 1 Introduction Multiscale simulation is a very important tool to reveal material behavior at different length and temporal scales through numerical simulations. Despite significant developments in materials simulation techniques, the goal of predicting reliably the properties and behaviors of new materials has not yet been achieved. This situation exists for several reasons that include a lack of full understanding of material behavior at different scales, absence of scaling laws, computational limitations, and difficulties associated with experimental measurements of material properties at micro and nano scales. For example, the information on the mechanical behavior of materials at nano level is not presently available as input to nanotechnology for the manufacturing of nanocomponents or microelectromechanical systems (MEMS). Scaling laws governing the mechanical behavior of materials from atomistic (nano), via mesoplastic (micro), to continuum (macro) scales are very important to numerous 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. Scaling laws are also important for 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 2 conventional indentation. Appropriate scaling laws may extend the extensive knowledge accumulated over time on material behavior at the macro (or continuum) level to the atomistic (or nano) level, via mesoplastic (micro) level. This work is focused primarily on developing coupling methodologies for continuum and atomistic simulations. At the continuum level, a new simulation method, the generalized interpolation material point (GIMP) method (Bardenhagen and Kober (2004)) is used. GIMP is developed based on the conventional material point method (MPM) (Sulsky et al. (1995)) with improved simulation stability. In order for GIMP to cover several length scales, and to couple MD simulations, refinement techniques must be developed. The thesis is divided into several chapters. Each chapter includes an introduction specific to that topic and a conclusion of findings. Chapter 2 describes the general theory of GIMP and multilevel refinement technique at the continuum level based on the Structured Adaptive Mesh Refinement Application Infrastructure (SAMRAI) (Hornung and Kohn (2002)). Parallel processing is implemented for GIMP based on SAMRAI and a contact algorithm is developed for the simulation of nanoindentation problems. Chapter 3 presents a structured mesh refinement algorithm for GIMP that is smoother than the method in Chapter 2. The new refinement technique is based on structured mesh and the refined mesh remains structured. The displacement boundary condition is explicitly introduced into the discretization of equation of motion in GIMP. Chapter 4 describes an application of the GIMP method for 3D simulations of bulk metallic glass (BMG) foams prepared by thermoplastic expansion. The material points are created based on the voxels in Xray tomography so that the microstructure of the foam is reconstructed exactly in GIMP. The simulation results of uniaxial compression 3 are compared with experimental ones. The densification of foam during compression is analyzed. Chapter 5 presents a coupling algorithm for GIMP and MD simulations. A new method to compute the atomic strain from MD simulation was introduced and validated. The coupling algorithm enables force, velocity, displacement, and energy compatibility at the handshaking region of the continuum and atomistic domains. Numerical simulations are performed to validate the coupling algorithm. The coupling algorithm Chapter 5 is extended to include discrete dislocations in Chapter 6. The discrete dislocation field is coupled with the continuum field using superposition of boundary conditions. The dislocations in the MD simulation can be detected and passed into the continuum region and vice versa. The summary and future works are presented in Chapter 7. The preliminary results of using GIMP to simulate metal cutting is presented in Chapter 7. This problem involves issues such as contact, large deformation, material failure, and adaptive mesh refinement. 4 Chapter 2 Multiscale Simulations Using Generalized Interpolation Material Point (GIMP) Method and SAMRAI Parallel Processing In the simulation of a wide range of mechanics problems including impact/contact/penetration and fracture, the material point method (MPM) (Sulsky, Zhou and Shreyer, (1995)) demonstrated its computational capabilities. To resolve alternating stress sign and instability problems associated with conventional MPM, Bardenhagen and Kober (2004) introduced recently the generalized interpolation material point (GIMP) method and implemented for onedimensional simulations. In this chapter GIMP has been extended to 2D and applied to simulate simple tension and indentation problems. For simulations spanning multiple length scales, based on the continuum mechanics approach, a parallel GIMP computational method is presented using the Structured Adaptive Mesh Refinement Application Infrastructure (SAMRAI). SAMRAI is used for multiprocessor distributed memory computations, as a platform for domain decomposition, and for multilevel refinement of the computational domain. Nested computational grid levels (with successive spatial and temporal refinements) are used in GIMP simulations to improve the computational accuracy and to reduce the overall computational time. The domain of each grid level is divided into multiple rectangular patches for parallel processing. This domain decomposition embedded in SAMRAI is very flexible when applied to GIMP. As an example to validate the parallel GIMP 5 computing scheme under SAMRAI parallel computing environment, numerical simulations with multiple length scales from nanometer to millimeter were conducted on a 2D nanoindentation problem. A contact algorithm in GIMP has also been developed for the treatment of contact pair between a rigid indenter and a deformable workpiece. GIMP results are compared with finite element results on indentation for validation. A GIMP nanoindentation problem with five levels of refinement was modeled using multiprocessors to demonstrate the potential capability of the parallel GIMP computation. 2.1 Introduction The material point method (MPM) has demonstrated its capabilities in addressing such problems as impact, upsetting, penetration, and contact (e.g. Sulsky, Zhou and Schreyer (1995); Sulsky and Schreyer (1996)). In MPM, two descriptions are used  one based on a collection of material points (Lagrangian) and the other based on a computational background grid (Eulerian), as proposed by Sulsky, Zhou and Schreyer (1995). A fixed structured mesh is generally used in the background throughout the MPM simulations. The material points are followed throughout the deformation of a solid to provide a Lagrangian description and the governing field equations are solved at the background grid nodes so that MPM is not subject to mesh entanglement. Compared to the finite element method (FEM), MPM takes advantage of both Eulerian and Lagrangian descriptions and possesses the capability of handling large deformations in a more natural manner so that mesh lockup problems present in FEM are avoided. Additionally, for problems involving contact, MPM provides a naturally nonslip contact algorithm to avoid the penetration between two bodies based on a common background mesh (Sulsky, Zhou and Schreyer (1995); Sulsky and Schreyer (1996)). One drawback of the 6 conventional MPM is that when the material points move across the cell boundaries during deformation, some numerical noise/errors can be generated, Bardenhagen and Kober (2004). To solve the instability problems associated with the conventional MPM simulations, Bardenhagen and Kober recently proposed the generalized interpolation material point method (GIMP) and implemented for onedimensional simulations. The present investigation extends the GIMP presented by BardenhagenKober to twodimensional simulations and applies it to simple tension and indentation problems. Furthermore, a refinement technique and a parallel processing scheme are developed so that the serial GIMP algorithm and code can be extended for parallel computation of large scale computations based on the continuum mechanics approach. Parallel processing has been used successfully in numerical analysis using different methods, such as FEM and boundary element method (BEM) (Mackerle (2003)) and molecular dynamics (MD) (Kalia and Nakano (1993)). The computational time on parallel processors can be reduced to a small fraction of the time consumed by a single processor at the same speed. Parallel processing generally involves issues, such as domain decomposition/partitioning, load balancing, parallel solver/algorithms, parallel mesh generation, and multigrid (Mackerle (2003)). Domain decomposition has been widely applied in parallel processing in FEM (Hsien (1997)). With partitioning of the overall computational domain, sequential FEM algorithm usually cannot be used directly in parallel processing without some modification, primarily due to the coupling of a large number of simultaneous linear equations. Remeshing is sometimes needed in each subdomain. The interfaces of neighboring subdomains must be meshed identically for subsequent communications (Mackerle (2003)). These problems are intrinsic to certain 7 numerical methods, such as FEM; however, they can be totally or partially avoided if other appropriate computational methods are used. For example, the domain decomposition is more straightforward for structured meshes, and large systems of coupled equations can be avoided, if explicit time integration is used. Recently a platform for parallel computation, namely, the structured adaptive mesh refinement application infrastructure (SAMRAI) (Hornung and Kohn (2002)), has been developed by the Center for Applied Scientific Computing at the Lawrence Livermore National Laboratory. SAMRAI has provided interfaces for userdefined data types so that material points carrying physical variables (mass, displacement, velocity, acceleration, stress, strain, etc.) can be readily defined. As a result, SAMRAI is very suitable for handling material points and their physical variables in MPM or its variant, GIMP. In this investigation SAMRAI is used for parallelizing GIMP. SAMRAI has also provided a foundation for parallel adaptive mesh refinement (AMR) with the use of either dynamic or static load balancing (Wissink, Hysom and Hornung (2003)). This function allows SAMRAI to process both spatial and temporal refinements in areas of interest, typically with high gradients in some physical variables (e.g., strains), and to use coarse mesh in the remaining areas. With the appropriate use of fine and coarse meshes in different regions, multiscale simulations using MPM can provide desired computational accuracy with reduced costs associated with computer memory and computational time. Material multiscale simulations span from electronic structure, atomistic scale, crystal scale, to macro/continuum scale (Horstemeyer, Baskes, Prantil, Philliber and Vonderheide (2003); Komanduri, Lu, Roy, Wang and Raff (2004)). Appropriate simulation algorithms can be used at various scales, e.g., ab initio computation for 8 electronic structure, molecular dynamics at the atomistic scale, crystal plasticity or mesoplasticity at the crystal scale, and continuum mechanics at the macro scale (Horstemeyer, Baskes, Prantil, Philliber and Vonderheide (2003)). At the continuum scale, FEM is generally used. Recently, the meshless local PetrovGalerkin (MLPG) method (Shen and Atluri (2004a, 2005)) and the continuum/lattice Green function method (Tewary and Read (2004)) have been used to couple with molecular dynamics seamlessly. The MLPG method is a simple, lesscostly alternative approach to FEM (Atluri and Shen (2002)). For the purposes of providing the insights into the discrete atomistic system and coupling with continuum, an equivalent continuum was defined in the MD region to compute the atomic stress based on the Smoothed Particle Hydrodynamics (SPH) method (Shen and Atluri (2004b)). The atomic stress tensor computed using the SPH method is more natural than other atomic stress formulations because it is in the nonvolumeaveraged form and rigorously satisfies the conservation of linear momentum. Hence, it is applicable to both homogeneous and inhomogeneous deformations. A tangent stiffness formulation was developed for both MLPG and MD regions and the displacements of the nodes and atoms are solved in one coupled set of linear equations. The MLPG/MD coupling has been demonstrated to be capable of enforcing the local balance equations in the handshaking region between continuum mechanics and molecular dynamics (Shen and Atluri (2005)). The simulation using parallel GIMP computing scheme in this investigation will focus on multiscales, e.g., from nanometer to millimeter, based on the continuum mechanics approach, namely, 2D GIMP. An example used for validating the simulation at several length scales at the continuum level is nanoindentation. It involves the contact issue 9 between a rigid indenter and a deformable workpiece. A contact algorithm, which allows the contact interface to be located in a few computational domains, is introduced in this study. The contact pressure is determined from solving a set of equations from multiple processors. Parallel GIMP results on nanoindentation are compared with FEM results using the ABAQUS/Explicit code. A nanoindentation model with five levels will be used; this model allows simulation from nanometer to millimeter scales. 2.2 Generalized Interpolation Material Point (GIMP) Method The governing equations in both conventional material point method (MPM) (Sulsky, Zhou and Sheryer (1995); Hu and Chen (2003); Bardenhagen (2002)) and generalized interpolation material point (GIMP) method (Bardenhagen and Kober (2004)) are briefly summarized in this section. The weak form of the momentum conservation equation in the conventional MPM is given by ∫ ⋅ Ω = −∫ ∇ Ω+ ∫ ⋅ + ∫ ⋅ Ω Ω Ω ∂Ω Ω ρw ad ρss : wd ρcs wdS ρw bs d , (21) where w is the test function, a is the acceleration, and ss, cs and b are the specific stress, specific traction, and specific body force, respectively. Ω is the current configuration and ∂Ω is the surface with applied traction. The material density, ρ, can be approximated as the sum of material point masses using a Dirac delta function Σ= = − Np p t p t M 1 ρ (x, ) δ (x x ) , where Np is the total number of material points and Mp is the mass of the material point. Upon discretization of Eq. (21) using the shape functions ( t ) i p N x , the governing equations at the background grid nodes become (see Sulsky, Zhou and Schreyer (1995)) ( )int ( t )ext i t i ti t i m a = f + f . (22) 10 where the lumped mass matrix is given by Σ= = N p p t p i p t i m M N 1 (x ) , (23) and the internal and external forces are given by Σ= = − ⋅∇ p t p N p i s t p p t i M N 1 ( )int , x f s , (24) Σ Σ = = = − − + p N p p t i p t p p N p t i p s t p p t ext i M h N M N 1 1 (f ) c , 1 (x ) b (x ) , (25) where h is the thickness of a boundary layer. At each time step, all variables for each material point, such as mass, velocity, and force are extrapolated to the grid nodes of the cell in which the material point resides. New nodal momenta are computed and used to update the physical variables carried by the material points. Thus, material points move relative to each other to represent deformation in a solid. A spatially fixed background grid is used throughout the MPM computation. MPM has already demonstrated its capabilities in solving a number of problems involving impact/contact/penetration. In case of large deformation, however, numerical noise, or errors have been observed, especially when material points have just crossed cell boundaries resulting in instability problems in the MPM simulations (see, e.g., Sulsky, Zhou and Schreyer (1995); Hu and Chen (2003); Bardenhagen and Kober (2004)). The primary cause for the problem has been attributed to the discontinuity of the gradient of the shape functions across the cell boundaries (see, e.g., Hu and Chen (2003); Bardenhagen and Kober (2004)). To resolve this problem, Bardenhagen and Kober (2004) proposed a generalized interpolation material point (GIMP) method. In GIMP, the interpolation between node i and material 11 point p is given by the volume averaged weighting function ∫ Ω∩Ω = p S d V S p i p ip (x) (x) x 1 χ , (26) where Vp is the current volume of the material point, (x) p χ is the characteristic function of the material point, and Si(x) is the node shape function. The role of the weighting function is the same as the shape function in conventional MPM. The modified equation of momentum conservation (Bardenhagen and Kober (2004)) can be written as Σ ∫ ∫ Σ ∫ Σ ∫ Ω∩Ω ∂Ω Ω∩Ω Ω∩Ω ⋅ + ⋅ ⋅ + ∇ = b v x c v x v x σ v x p d d V m d d V p p p p p p p p p p p p p p δ δ χ δ χ δ χ : & , (27) where δv is an admissible velocity field, p p& is the rate of change of the material point momentum. Eq. (27) can be further discretized and solved at the grid nodes, Bardenhagen and Kober (2004). Herein, the weighting function S ip is C1 continuous under the spatially fixed background grid. Consequently the noises associated with material point crossing cell boundaries in the conventional MPM can be minimized. In this chapter, the GIMP presented by BardenhagenKober is implemented for twodimensional simulations. a refinement technique and a parallel processing scheme to extend the serial GIMP algorithm to code large scale parallel computing have been developed. The capability of parallel GIMP computing has been demonstrated by modeling nanoindentation problem. A contact algorithm has been developed to address the contact problem between the rigid indenter and the deformable workpiece. Next, the contact algorithm developed in this investigation is described. 12 2.2.1 Contact algorithm in GIMP Nanoindentation involves a contact pair of a rigid indenter and a deformable workpiece. The contact interaction between these two surfaces is governed by the Newton’s third law and Coulomb’s friction law as well as the boundary compatibility condition at the contact interface (Oden and Pires (1983); Zhong (1993)). While MPM can prevent the penetration at the interface automatically, it uses a single mesh for the two bodies. At the contact surface, all components of the variables are interpolated to the nodes from both bodies using Eqs. (23)(25). As a result, MPM using a single mesh tends to induce early contact in approaching and late separation when two parts move away from each other. So, MPM cannot model the contact behavior between two parts correctly. Hu and Chen (2003) proposed a multimesh MPM algorithm to release the noslip constraint inherent in the MPM using a single mesh. In the multimesh MPM, in addition to a common mesh for all objects, there is an individual mesh for each of the objects under consideration. All meshes are identical, i.e. nodal locations are the same. The multimesh can be generated by creating multiple nodal fields for each node. Each nodal field corresponds to an object. In multimesh MPM scheme, the nodal masses and forces are mapped from the material points of each object to its own mesh. The nodal values are transferred to the corresponding nodes in the common mesh. When the values at a node of the common background mesh involve contributions from two parts, the contact between two parts occurs so that this node is defined as an overlapped node. Otherwise, two parts move independently. This multimesh algorithm can handle sliding and separation for the contact pair. However, in using the multimesh for contact problem in GIMP, the interaction at the overlapped nodes is still activated too early before the actual contact of 13 the material points occurs. Fig. 21 is an example illustrating early contact when Part 1 is moving toward Part 2. The four bottom particles of Part 1, labeled in hollow circles, have come into the cells of Part 2, the nodes of which cells are labeled in three dashed circles. Physical variables (e.g., normal force, and velocity component normal to the contact surface) in Part 1 will be interpolated onto these three overlapped nodes. The physical variables in the three overlapped nodes will be further interpolated into material points within the top layer of cells in Part 2, and contribute to the stress and deformation in the entire Part 2. With this treatment in the previous multimesh algorithm, even though Parts 1 and 2 are not in physical contact, the particles in Part 2 will contribute to the physical variables of particles in Part 1, leading to numerical early contact, and vice versa, through the overlapped nodes. Similar situation occurs when Part 1 is retracting from Part 2, resulting in late separation of two parts. Unless other measures are taken to prevent these physically incorrect early contact and late separation problems, they could cause large errors in GIMP and must be corrected in contact problems. Fig. 21: Illustration of early contact in multimesh In this chapter, a new contact algorithm is developed for GIMP simulations. Fig. 22 illustrates the contact algorithm for the contact pair between a rigid indenter and a Part 1 Part 2 Overlapped nodes 14 deformable workpiece. Although circular points are used in this schematic diagram, it should be noted that the points are representations of areas occupied by these points, based on the GIMP algorithm. A frictionless contact is assumed in this investigation. At the beginning of a time step, a material point is located at point A. At the end of this time step, the material point moves to B, if there is no contact interaction. Fig. 22: Schematic of contact algorithm between the rigid indenter and the deformable workpiece To satisfy the displacement compatibility condition, the material point has to be brought to the indenter surface and kept in contact with the indenter. The contact velocity correction c p V can be determined based on the rigid surface orientation indicated by its unit outward normal vector n. The final location of the material point is set to C by a contact pressure. Hence, the velocity of a material point p under contact can be determined by Σ Σ Σ Σ Δ + + Δ = + + Δ = = N i ip i c i N i ip i i i N i ip i c i i i N i p i ip S m S t m t S m S t p F F V V p F F 0 0 0 ( 0 ) , (28) where N is the number of nodes contributing to material point p, c i F is the contact force on node i, 0 i p and 0 i F are the nodal momentum and force without consideration of the A C B n D 0 p V c p V p V Rigid surface Rigid indenter Indenter side Workpiece side Workpiece 15 contact, respectively. The velocity 0 p V of the material point without the consideration of contact is given by Σ + Δ = N i ip i i i p S m 0 0 t V0 p F . (29) The contact force c i F on node i is the resultant of the contact pressures on the neighboring particles, and can be computed in terms of contact pressures using the approach given by Bardenhagen and Kober (2004), i.e., Σ ∫ = ∂Ω = Q q c i q c i q S ds 1 F (x)P , (210) where Q is the total number of material points in contact with the indenter. If the contact pressure c q P is assumed to be constant in the contact area occupied by material point q, then, Σ= = Q q c iq q c i T 1 F P , where ∫ ∂Ω = q T S ds iq i (x) . Since c p p p V = V0 + V , the contact velocity c p V for material point p is given by = Δ Σ Σ N i Q q c iq q i c ip p T m V t S P . (211) Eq. (211) can be established for each material point in contact. At each material point there is an unknown contact pressure c q P . Therefore, the number of unknown pressures, c q P , is equal to the number of points in contact. In parallel computing, points in contact might be located in different domains processed by different processors. Consequently, a parallel solver is needed to solve Eq. (211) in this investigation. Since contact can only occur on the outer surface of an object, Eq. (211) is solved analytically under the physical contact condition Pc ⋅n < 0 q to find the contact pressure at all material points in 16 contact with the indenter. The contact pressure is then extrapolated to the nodes from the contact material points to update the total nodal forces. 2.2.2 Numerical implementation Fig. 23: Material points in cells and the weighting function in 2D GIMP Considering the case where initially there are four material points in a cell, the 2D weighting function is depicted in Fig. 23. To compute the weighting function, χp( x ) is one in the current region occupied by the material point p and zero elsewhere. In this figure, one node is at the origin and the horizontal axes give material point positions normalized by the cell size. Fig. 23 is based on the same material point characteristic function and node shape function as in Bardenhagen and Kober (2004). It is noted that the computation of the weighting function in the deformed state involves some practical difficulties because the integration boundaries in Eq. (27) can be difficult to obtain. To circumvent this problem, it is assumed that the shape of the region occupied by the four material points remains rectangular without rotation, so that Eq. (26) can be evaluated analytically. This assumption leads to significant saving in the computational time while introducing only small errors. Using this assumption, GIMP is extended to 2D simulations and the results are presented in Section 4. 1 0.5 0 0.5 1 1.5 1.5 1 0.5 0 0.5 1 1.5 0 0.2 0.4 0.6 0.8 17 2.3 Parallel Computing Scheme Using GIMP with SAMRAI 2.3.1 SAMRAI The Structured Adaptive Mesh Refinement Application Infrastructure (SAMRAI), a scientific computational package for structured adaptive mesh refinement and parallel computation, is used with the GIMP for parallel computation of largescale simulations. SAMRAI is chosen because of its similarity in grid structure with GIMP. In GIMP, the computation is usually independent of the background grid mesh so that a structured spatiallyfixed mesh can be used throughout the entire simulation process. This advantage makes GIMP highly suitable for parallel computation, as the domain decomposition for structured mesh can be easily performed and no remeshing is required. Thus the complexity and inefficiency associated with parallel processing can be avoided. In SAMRAI, the computational domain is defined as a hierarchy of nested grid levels of mesh refinement (Berger and Oliger (1984)) as shown in Fig. 24. Each grid level is divided into nonoverlapping, logicallyrectangular patches, each of which is a cluster of computational cells. Indices are used extensively in SAMRAI to manage grid levels and patches. For example, patch connectivity is managed by the cell indices. The organization of the computational mesh into a hierarchy of levels of patches allows data communication and computation to be expressed in geometricallyintuitive box calculus operations. Communication patterns for data dependencies among patches can be computed in parallel without interprocessor communications, since the mesh configuration is replicated readily across processor memories. Interprocessor communications, i.e., data communications between patches on the same as well as neighboring levels, are predefined by SAMRAI communication schedules. Problem 18 specific communication interfaces are also provided by SAMRAI. SAMRAI supports several data types defined in a patch, such as cellcentered data, nodecentered data, and facecentered data. These data are stored as arrays to allow numerical subroutines to be separated easily from the implementation of mesh data structures. Userdefined data structures over a patch, which can be accessed through cell index, are supported by SAMRAI. These characteristics make SAMRAI a very flexible parallel computing environment for numerous physics applications (Wissink, Hysom, and Hornung (2003)). Fig. 24: Illustration of a hierarchy of three nested grid levels of mesh refinement 2.3.2 Spatial and temporal refinements In the application of SAMRAI to largescale GIMP simulations, the techniques for refinement, both spatial and temporal, have to be developed to achieve high accuracy in areas of high stress/strain gradients while reducing the overall computational time by using coarse mesh in regions of low stress/strain gradients. Since a structured mesh is used in GIMP, the refinement can be implemented by imposing fine levels of subgrids at locations of interest, using the approach adopted by Berger and Oliger (1984) in SAMRAI. The scheme for the structured grid refinement is illustrated in Fig. 24. The Level 1 Level 2 Level 3 19 cell size ratio, also called the refinement ratio, of two neighboring levels is always an integer for convenience. The advantage of this refinement technique is that nesting relationships between different levels can be handled. A material point in GIMP can be split into several small material points. Tan and Nairn (2004) proposed a criterion to split material points based on local deformation gradient. If the refinement ratio is two in each direction, one coarse material point can be split into four material points in the next fine level in 2D case. However, this splitting technique can become complicated when conservation of energy and momentum have to be enforced. In this chapter, a more natural refinement approach is developed to avoid direct splitting and merging processes by using material points of the same size and mass in overlapped region (called ghost region) between two neighboring levels. Fig. 25 shows two neighboring coarse and fine grid levels in 2D GIMP computations with a refinement ratio of two. The thick line represents the physical boundary of the fine level with four layers of ghost cells. Initially, four material points are assigned to each cell at the fine level. At the coarse level, the portion overlapped by the fine level is assigned with 16 material points per cell. Hence, these material points have the same size and initial positions as those at the fine level. The rest of the coarse level is assigned with four material points per cell. GIMP provides a natural coupling of the material points with different sizes at the same grid level. This is because the weighting function depends on the characteristic size of the material points and cell length and the interpolation between the nodes and material points is weighed by the mass of the material point. In GIMP computation, each level is computed independently with the physical variables communicated through the ghost regions between neighboring levels. Two data exchange 20 processes, namely, refinement and coarsening are used in the communication. Refinement process passes information from the coarse level to the immediate fine level, while coarsening process will pass information from fine level to the next coarse level. In the refinement process, physical variables at the fine material points inside the thick lines in Fig. 25 are copied directly to replace the material points at the coarse level. In the coarsening process, the physical variables at coarse material points are copied to the ghost cells of the immediate fine level. In the refinement, the material points located in the ghost cells at the fine level are eliminated first, and the material points in the corresponding region of the coarse level (with the same size as points in the immediate fine level) are copied to ghost cells at the fine level. In the copying process, if some (small) material points in the coarse level fall into the interior (inside the square on the right of Fig. 25) of the immediate fine level, these points will not be copied to this fine level, as points in the fine level will be able to carry over all computations in the interior of the fine level already. In coarsening, the material points of the coarse level located in the region overlapped by the fine level (inside the square on the left of Fig. 25) are eliminated first, and the material points in the fine level are then copied to the immediate coarse level. With these refinement/coarsening operations, the material points can move around freely, including moving outside the original level in large deformations. To ensure that this coarsening process can still be performed reliably during deformation, sufficiently wide region of cells should be assigned with refined material points at the coarse level so that the ghost cells of the fine level always stay within the region with fine material points on the coarse level. At the coarse level, the interior cells covered by the fine level do not participate in 21 the computation and there are no material points inside (see Fig. 25). Fig. 25: Two neighboring coarse and fine grid levels in 2D GIMP computations Fig. 26: Nested multilevel refinement and reduction in the number of cells with number of levels The refinement techniques can be applied for multiple times at the regions of interest, such as the stress concentration regions. A fixed refinement ratio of two between two neighboring levels is very effective in reducing the total number of computational cells. Fig. 5 shows nested multilevel refinement and its corresponding relation between the total number of cells and the number of grid levels. The cell percentage represents the ratio of the total number of cells with multilevel refinement mesh to the total number of cells with onelevel finest mesh. If each fine level occupies one quarter of the neighboring coarser level, as shown in Fig. 26 (a), the cell percentage as a function of Coarse Fine Level 1 Level 2 3 (b) 0 20 40 60 80 100 1 2 3 4 5 6 7 Number of levels Cell percentage (%) (a) 22 the number of grid levels can be calculated, as shown in Fig. 26 (b). For example, when totally four levels of successive refinements are used the total number of cells is about 8% of that of one uniform fine mesh. A reduction in the number of computational cells leads to a reduction in the number of material points. Hence, the total amount of computational time can be reduced significantly. However, refinement and coarsening communications will cost additional computational time, as will be discussed in Section 4. Another advantage of the multilevel refinement is that it allows for temporal refinement. Since the computation at each grid level is conducted independently, different time step increments can be used for computation at different levels. For example, a smaller time step increment can be used for the fine level to improve computational accuracy, while a larger time step increment can be used for the coarse level. Since the refinement ratio is an integer, the time step increment ratio should also be an integer for convenience in the computation and data communication/synchronization. For example, in Fig. 25, when the refinement ratios in both directions are fixed at two, the time step increment ratio should be set to two as well. As a result, two time step computations are performed at the fine level, and results are passed over to the immediate coarse level to couple with the results at the coarse level. 2.3.3 Domain decomposition GIMP uses structured mesh, consistent with SAMRAI, so that domain decomposition is straightforward and no remeshing, in general, is necessary. Fig. 27 (a) shows a twodimensional computational domain decomposed into two patches separated by a horizontal dash line. The elliptical solid object with different boundary conditions applied at different regions is inside this domain/grid. After discretization, there are a certain 23 number of material points and part of the boundary in a patch, which will be computed individually. It may be noted that patch boundary does not have to coincide with the boundary of the material continuum. The patch boundary is always chosen to be larger than the region occupied by the material continuum so that there is extra space for the material to deform. This will not cause any additional computational burden as the GIMP computation is only carried out on material points inside the patch. Each patch can be processed by a single processor and the convenience in creating patches will provide great flexibility in parallel processing. Fig. 27: A computational domain of two patches in one grid level Communication between two neighboring patches is realized through information sharing in the region overlapped by the two patches. The overlapped regions are also called “ghost” regions, as shown in Fig. 27 (b). The ghost cells are denoted by dash lines. For ease of visualization, only the ghost cells overlapped by the other patch are shown and the ghost cells along the other three sides of a patch are not shown. On one grid level, patches can communicate with each other by simply copying data from one patch to another at the same computation time (Fig. 27 (b)). Using the material point information from the previous time step, and the physical boundary conditions, each patch is ready to advance one more time step. At this time, the material point information in the outermost Patch 2 Patch 1 Ghost of patch 1 Copy from patch 2 Ghost of patch 2 Copy from patch 1 (a) (b) 24 layer of the ghost cells becomes inaccurate. For instance, one outermost grid node in patch one, marked by the circle, obtains information from eight material points before advancing to the next step. After advancing, it extrapolates to eight material points. However, in patch two, the grid node at the same location obtains information from sixteen surrounding material points. It extrapolates to these sixteen material points after advancing. Typically, after each step, the material points in the next inner layer in the ghost region become inaccurate as well. Ghost cells and material points are attached to each patch to ensure accuracy of the interior. Each patch can be computed independently for one GIMP step since the momentum conservation equation is solved at each node and there are no coupled equations to solve. No data exchange is necessary during the GIMP step. Therefore, different patches can be assigned to different processors for parallel processing. After one GIMP step, the data in the ghost cells will be updated. Fig. 28: Flowchart showing advancement of grid levels recursively starting from the coarsest to the finest level in GIMP AdvanceLevel (level_number) Advance each patch in level by level time step Start level_number=0 Next finer level exists AND level time > next finer level time? level_number = level_number+1 Data communication for level YES NO End 25 Copying material points to ghost cells involves data exchange between processors, which costs additional time. The more the number layers of ghost cells, the longer the time needed for communication, but communication can be performed less frequently. A minimum of two layers of ghost cells are necessary to ensure that computation at the material points inside a patch is always correct. If three levels of ghost cells are chosen, the communication can be performed after every two increments of each patch. With these refinements and domain decomposition schemes for GIMP, it is possible to implement GIMP into the SAMRAI platform. In this study, the refinement ratio is chosen as two. Four layers of ghost cells are augmented to a patch such that data communications, including both data exchange on the same level and between neighboring levels, are performed every two timesteps for each fine grid level. This is critical because data exchange between levels has to be performed when the two levels are synchronized. Fig. 28 shows the flowchart advancing all grid levels recursively starting from the coarsest level for one coarsest time step. It may be noted that the sequential GIMP algorithm can be used to advance each patch without modification. 2.4 Numerical Examples and Results A Beowulf Linux cluster of 8 identical PCs were used in the simulations. Each PC has a Pentium 4 processor with a 2.4 GHz CPU, 512 MB RAM except that the master node has a memory of 1 GB. A gigabit switch is used to connect the network. Two examples are used for validation of the 2D parallel GIMP computing under SAMRAI platform. The first example is simple tension of polycrystalline silicon under plane strain conditions. The material is assumed to be homogeneous, isotropic, and linear elastic. The Young’s modulus is 170 GPa and the Poisson’s ratio is 0.18. One end is 26 constrained along the X direction while a normal traction is ramped up on the other end. The size of the tensile model is 0.06 mm × 0.04 mm. The length of a square grid cell is 0.002 mm and the time step is 5×104 ps. For verification, the same problem was simulated using both conventional MPM and FEM (ABAQUS/Explicit). Fig. 29 shows GIMP, MPM and FEM simulation results of normal stresses in the Xdirection at different increments from a simple tension problem. The simulation using the conventional MPM in Fig. 29 (a) shows material separation close to the free end with severe numerical instability after 275 increments. Fig. 29 (b) and (c) show the normal stress distribution in the tensile direction and deformation after 500 increments from FEM and GIMP. It may be noted that FEM results show a stress contour plot on the deformed mesh while the GIMP results show a discrete scattered plot of material points. These two results are in good agreement with the difference in the maximum value being less than 10%. Fig. 29: Simulation results of tensile stress contours for a simple tension problem (c) FEM (a) Conventional MPM (b) GIMP 27 Fig. 210: Loading conditions for a simple indentation problem In the second example, indentation on the same silicon material is simulated. The workpiece is subjected to a pressure applied in the middle of the top surface (Fig. 210 (a)) under plain strain conditions with a thickness of 0.001 mm for computing the mass and forces. The magnitude of the pressure increases linearly with time for the first 1500 increments, and is then kept constant (see Fig. 210 (b)). The cell size is 0.001 mm in both directions and the time step is 20 ps for both FEM and GIMP simulations. Due to symmetry, only one half of the workpiece is modeled. This simulation is performed with two patches in one uniform grid level. Two processors are used and one patch is assigned to each processor. Fig. 211 shows GIMP and FEM results of normal stresses in the Ydirection at different increments. The dashed line in Fig. 211 (a) is the boundary between the two patches. Fig. 211 (a) and (b) are plots of normal stresses in Ydirection at 500 time increments for GIMP and FEM simulations. The difference in stress values in Fig. 211 (a) and (b) is less than 5%. It should be noted that the FEM simulation aborted at 1348 increments due to excessive element distortion. The GIMP simulation did not encounter this problem. Fig. 211 (c) shows the GIMP stress result after 2000 increments. This demonstrates the capability of GIMP in handling excessive distortions. X Y (a) (b) Time Pressure t = 3×10−8 s 60 GPa 28 Fig. 211: GIMP and FEM results of normal stress variation in the Ydirection at different increments In order to validate the multilevel refinement algorithm and parallel communication, as well as the proposed contact algorithm, a simulation of nanoindentation with a wedge indenter was conducted under 2D plane strain conditions. The workpiece is aluminum and the indenter is assumed to be rigid. Fig. 11 shows the indentation model. The area below the indenter where high stress gradients are expected is refined, as shown in Fig. 212 (a). A prescribed velocity was applied on the indenter, as shown in Fig. 212 (b). (c) GIMP at 2000 increments (a) GIMP at 500 increments (b) FEM at 500 increments 29 The work piece dimensions are 60 μm × 40 μm. It is fixed in the Ydirection at the bottom. Only half of the model is simulated because of symmetry. The cell sizes are 500 nm, 250 nm and 125 nm for levels 1, 2 and 3, respectively. Each level is divided into four patches with approximately the same size. The maximum indentation depth in the simulation is about 450 nm. The dotted lines in Fig. 212 (a) illustrate the four patches in level 1. For comparison, an explicit FEM simulation (using ABAQUS/Explicit) was carried out under the same conditions. The FEM element size is uniform and is the same size as the finest GIMP background grid size. In this example, the maximum indentation depth (450 nm) is relatively small compared to the finest cell size, so that FEM simulation has not encountered excessive mesh distortion. Fig. 212: Schematic of 2D indentation showing (a) three levels of refinement and (b) the indenter velocity history Fig. 213 shows a comparison of contours of normal stresses in the Ydirection at the maximum depth for both FEM and parallel GIMP simulations. The axis of symmetry of the workpiece is located at X = 0.03 mm. For FEM, the plot is the contour of nodal stresses with deformed positions, and for GIMP, it is a discrete scattered plot of stress at deformed material points. The area below the indenter with high stress gradients is refined as shown in Fig. 212 (a) for parallel GIMP computation. The borders of grid levels 2 and 3 can barely be seen in Fig. 213 (b) due to the use of high material point (a) (b) Time V 20 m/s V Level 1 Level 2 3 X Y 30 density. Fig. 214 is a closeup view of shear stresses in which the three grid levels are shown. Results show that the normal and shear stresses from both parallel GIMP and FEM simulations agree very well. The difference of the normal stresses in the Ydirection for the material point in contact with the indenter tip and the stress of the FEM node at the same location is 4.4%. It may be noted that some nonsmoothness in the GIMP stresses around the level boarders can be seen. This nonsmoothness is caused by the refinement and coarsening and the error associated with this is negligible for these simulations. Fig. 213: Comparison of normal stresses in Ydirection from FEM and GIMP GIMP simulations using a uniform cell size of 500 nm and 125 nm were performed under the same conditions as in Fig. 212 to further verify the refinement/coarsening algorithm. Fig. 215 shows normal stresses in the Ydirection and shear stresses in GIMP simulations using 500 nm uniform cells. In this case the material points in contact with the indenter are 16 times (4 times in each direction) larger than those with two levels of refinements. In general, the stress magnitudes agree with those in Fig. 213 and Fig. 214. Fig. 216 (a) shows indentation load versus depth curves from FEM and GIMP with (a) FEM (b) GIMP 31 different grid sizes. From GIMP computations, the load versus depth curves with three levels of refinement agree very well with the results from a uniform finest mesh. The load versus depth curve from the FEM simulation with a uniform cell size of 125 nm under the same boundary conditions is plotted for comparison. It can be seen from Fig. 216 (a) that the trend of the load versus depth plots from FEM and GIMP simulations are similar. The difference in indentation load at the end of loading, which corresponds to 450 nm of indentation depth, is 5.9% between FEM and GIMP with 3 grid levels. When the depth is less than 100 nm, there is only one material point in contact with the rigid indenter. The assumption of constant pressure causes large differences under this circumstance. However, if the GIMP cell size is further refined to 62.5 nm and the size of the material point is 31.25 nm, the difference between GIMP and FEM becomes smaller, as can be seen in Fig. 216 (b). Fig. 214: Comparison of shear stresses from FEM and GIMP Other simulations were conducted for the same problem with three levels of refinement using different number of processors to test the efficiency of parallel computing. The number of patches at each level is the same as the number of processors and the size of (a) FEM (b) GIMP 32 each patch is approximately the same. The resultant stress distribution and indentation load versus depth plots are the same as the previous results. The average time per computational step is 7.14 s when one processor is used and is reduced to 4.26, 3.40, 2.18 s, respectively when two, three, and four processors are used. When four processors are used, the CPU time per step is only 30.5% of that of one processor. This gives a speedup by a factor of 3.28. In the ideal case without communication overhead, the speedup would be 4. The reduction in speedup from the ideal number is because of the time involved in data communication between processors. It has been observed that the refinement and coarsening algorithm consume most of the communication time. Moreover, in refinement and coarsening, most of the time is taken to search for the corresponding material point in another grid level. This portion of the computational time can be reduced, if improved searching algorithm or more optimized algorithm for the storage of material points can be implemented. Fig. 215: Normal and shear stresses of GIMP simulations with a uniform cell size of 500 nm The manual refinement for the indentation problem is adequate since the region of high (a) Normal Stress (b) Shear Stress 33 stress gradient is known to occur below the indenter. The finest level covers the indenter and part of the specimen. With the same initial condition, the results at the finest level is identical to the results in the same area if a uniform fine mesh is used for the entire domain that requires much longer computational time. The computational load of each processor is balanced statically by assigning approximately the same number of material points to each processor. Dynamic load balance is supported by SAMRAI and can potentially improve the efficiency of the simulation. Fig. 216: Indentation load versus depth curves from FEM and GIMP with different grid sizes To demonstrate the capability of the algorithm developed in this investigation for multiscale simulation, an indentation model with multiple length scales is simulated with eight processors. The dimensions of the workpiece are 0.25 mm × 0.125 mm. Initially, the velocity of the indenter increases from 0 to 150 m/s linearly with time and is then kept constant. Five successive levels of refinement are used in this simulation. The smallest material point represents an area of 64 nm × 64 nm, and the largest material point covers an area of 1 μm × 1 μm. Each level is divided into 8 equalsized patches for best load balance. Since the contact surface can evolve into several patches, a parallel 0 1 2 3 4 5 6 0 100 200 300 400 500 Depth (nm) Load (mN) GIMPMcoarse (500 nm) GIMPMfine (125 nm) GIMPM3 levels FEMfine (125 nm) 0 0.5 1 1.5 2 2.5 0 40 80 120 160 Depth (nm) Load (mN) GIMPM500 nm GIMPM125 nm GIMPM62.5 nm FEM62.5 nm (a) (b) 34 solver is implemented to solve Eq. (211) to find the contact pressure based on the Portable, Extensible Toolkit for Scientific Computation (PETSc). An aluminum workpiece is chosen with the Young’s modulus and Poisson’s ratio of 70 GPa and 0.33, respectively. The maximum indentation depth was 9.8 μm in this simulation (i.e., 153 times the size of the finest material point). It took nine hours to simulate this problem with eight processors. Fig. 217 (a) gives the normal stress distributions and Fig. 217 (b) shows normal stress distribution for the finest two levels. The relative large deformation in this multiscale nanoindentation problem could not be handled by FEM due to excessive distortion in the FEM mesh. However, the parallel GIMP code was able to complete the entire loading/unloading processes without any difficulty. This example shows clearly the advantage of GIMP for multiscale simulations over FEM. Fig. 217: Multiscale simulation of nanoindentation with five levels of refinement 2.5 Conclusions The following are specific conclusions based on the results of this investigation: A 2D generalized interpolation material point (GIMP) method has been implemented to address problems, such as particle flyingoff and alternating stress sign associated with (a) (b) 35 conventional MPM in case of relatively large deformation. To conduct multiple length scale simulations, a parallel computing scheme has been presented using GIMP under SAMRAI parallel computing environment in which multilevel grids are used for spatial and temporal refinements. A refinement/coarsening algorithm, based on material points of GIMP in two grid levels, has been developed for communication between neighboring grid levels of different refinements. With increase in the refinement levels, as well as decrease in the time step increments, the computational accuracy is greatly improved in the region of interest while the overall computational time is reduced. The computation at each grid level is performed recursively to ensure that the refinement and coarsening are performed when the two neighboring levels are synchronized. 2D MPM and GIMP were applied to simple tension and indentation problems to validate the GIMP algorithm. GIMP results agree very well with FEM results for these two examples provided that the deformations are small. The noise and instability problems present in conventional MPM are not observed in the GIMP simulations. As the deformation is increased, GIMP continued to execute while FEM aborted due to element distortion. Also GIMP results are stable. Thus GIMP is able to handle relatively large deformation problems. For the nanoindentation problem, a GIMP algorithm for the contact between a rigid indenter and a deformable workpiece was developed. A reasonably good agreement between GIMP and FEM results was reached, validating the contact algorithm presented in this investigation. Another nanoindentation example with multiple length scales from a few nanometers to 36 submillimeters was simulated and numerical results validated the parallel GIMP computing with the use of SAMRAI. 37 Chapter 3 Structured Mesh Refinement in GIMP Method for Simulation of Dynamic Problems The generalized interpolation material point (GIMP) method, recently developed using a C1 continuous weighting function, has solved the numerical noise problem associated with material points just crossing the cell borders, so that it is suitable for simulation of relatively large deformation problems. However, this method typically uses a uniform mesh in computation when one level of material points is used, thus limiting its effectiveness in dealing with structures involving areas of high stress gradients. In this chapter, a spatial refinement scheme of the structured grid for GIMP is presented for simulations with highly localized stress gradients. A uniform structured background grid is used in each refinement zone for interpolation in GIMP for ease of generating and duplicating structured grid in parallel processing. The concept of influence zone for the background node and transitional node is introduced for the mesh size transition. The grid shape function for the transitional node is modified accordingly, whereas the computation of the weighting function in GIMP remains the same. Two other issues are also addressed to improve the GIMP method. The displacement boundary conditions are introduced into the discretization of the momentum conservation equation in GIMP, and a method is implemented to track the deformation of the material particles by tracking the position of the particle corners to resolve the problem of artificial separation of material particles in GIMP simulations. Numerical 38 simulations of several problems, such as tension, indentation, stress concentration and stress distribution near a crack (mode I crack problem) are presented to validate this refinement scheme. 3.1 Introduction The material point method (MPM) uses a collection of material points, mathematically represented by Dirac delta functions to represent a material continuum (Sulsky, Zhou, and Schreyer (1995); Hu and Chen (2003); Guilkey and Weiss (2003)). A spatially fixed background grid, and interpolation between grid nodes and material points are introduced to track physical variables carried by the material points in the Lagrangian description. Field equations are solved on the background grid in the Eulerian description. Physical variables are interpolated from the solutions on the background grids to material points back and forth for solution and convection of physical variables. In general, the isoparametric shape functions, same as those used in the finite element method (FEM), are used. As the MPM simulation is independent of the background grid, a structured grid is usually employed for purposes of simplicity. The movement of the material points represents the deformation of the continuum. MPM has been demonstrated to be capable of handling large deformations in a natural way (Sulsky, Zhou, and Schreyer (1995)). However, primarily due to the discontinuity of the gradient of the interpolation function at the borders of the neighboring cells, artificial noise can be introduced when the material points move just across the grid cell boundaries, leading to simulation instability for MPM. The generalized interpolation material point (GIMP) method, introduced by Bardenhagen and Kober (2004) can resolve this problem. In GIMP a C1 continuous interpolation function is used and each material point/particle occupies a region. GIMP 39 has been demonstrated to be stable and capable of handling relatively large deformations (Ma, Lu, Wang, Roy, Hornung, Wissink, and Komanduri (2005)). The current MPM typically uses a uniform background mesh for solving the field equations. However, this is not efficient when stress gradients are high such as stress concentrations in a plate with a hole, or the stress field of a workpiece under indentation. In contrast, transitional mesh is effective in solving problems involving rapidly varying stress in an area. Wang, Karuppiah, Lu, Roy, and Komanduri (2005) have presented a method using an irregular background mesh to deal with problems involving rapidly varying stress, such as stress field near a crack. However, this approach does not use regular structured background mesh so that mesh generation encounters the same difficulty as FEM, and leads to the loss of some advantage of MPM on the ease of generating mesh for a complex problem. The use of structured grid in GIMP has facilitated the implementation of GIMP in parallel processing. A refinement scheme based on splitting and merging material particles was proposed by Tan and Nairn (2002). Recently, a multilevel refinement algorithm has been developed for parallel processing using the structured adaptive mesh refinement application infrastructure (SAMRAI) (Hornung and Kohn (2002); Ma, Lu, Wang, Roy, Hornung, Wissink and Komanduri (2005)). The computational domain is divided into multiple nested levels of refinement. Each grid level is uniform but has a different cell size. Smaller material particles and smaller cell sizes are used in each finer level. Two neighboring levels are connected by overlapped material particles of the same size and data communication between levels is performed at predefined intervals. However, the refinement through material particles requires extra communication and 40 simulation time. In this chapter, a refinement for GIMP based on the transitional grid nodes is developed. This refinement is natural and does not involve extra simulation time. Moreover, the refined grid remains uniformly structured in each refinement region. While the problem associated with artificial noise has been resolved with the use of GIMP method, it has been observed recently that material separation could occur if the deformation of the material particles was not tracked (Guilkey (2005)). Tracking the deformation of material particles properly in GIMP is necessary especially when the material particles are stretched. In this chapter, an approach is developed for tracking the particle deformation to resolve material point separation problem. This chapter focuses on the refinement scheme for structured grid. Several numerical problems, such as tension, indentation, stress concentration and stress distribution near a crack (mode I crack problem) were simulated to verify this refinement algorithm, as well as to demonstrate the effectiveness of tracking particle deformations. 3.2 GIMP For the purpose of completeness, the basic equations in GIMP (Bardenhagen and Kober (2004)) are summarized here. In dynamic simulations, the mass and momentum conservation equations are given by + ρ∇⋅ v = 0 ρ dt d , and (31) ρa = ∇⋅σ + b in Ω, (32) where ρ is the material density, a is the acceleration, σ and b are the Cauchy stress and body force density, respectively. The displacement and traction boundary conditions are given as 41 u = u on u ∂Ω , (33) τ = τ on τ ∂Ω , (34) where ∂Ω ⊂ ∂Ω, u ∂Ω ⊂ ∂Ω τ and ∂Ω ∩∂Ω = 0. u τ In variational form, the momentum conservation equation can be written as ∫ ∫ ∫ ∫ Ω Ω Ω ∂Ω ⋅ = ∇ ⋅ + ⋅ − − ⋅ u ρa δvdx σ δvdx b δvdx α (u u) δvdx , (35) whereδv is an admissible velocity field, α is a penalty parameter introduced herein to impose the essential boundary conditions and α >>1 (Atluri and Zhu (1998); Atluri and Zhu (2000)). Applying the chain rule, (∇⋅σ) ⋅δv = ∇⋅ (σδv) − σ :∇δv , and the divergence theorem, Eq. (35) can be written as ∫ ∫ ∫ ∫ ∫ ∫ ∂Ω ∂Ω Ω Ω Ω ∂Ω + ⋅ − − ⋅ ⋅ + ∇ = ⋅ + ⋅ u u dS dS d d d dS u τ v u u v a v x σ v x b v x τ v δ α δ ρ δ δ δ δ τ ( ) : , (36) where u τ is the resultant traction due to the displacement boundary condition on u ∂Ω . In GIMP, the domain Ω is discretized into a collection of material particles, with p Ω as the domain of particle p. The physical quantities, such as the mass, stress and momentum can be defined for each particle. For example, the momentum for particle p can be expressed as ∫ Ω = p d p p p ρ (x)v(x)χ (x) x , where v( x ) is the velocity and ( ) p χ x is the particle characteristic function. The momentum conservation equation can be discretized as Σ ∫ Σ ∫ Σ ∫ Σ ∫ Σ ∫ Σ ∫ ∂Ω ∩Ω ∂Ω ∩Ω ∂Ω ∩Ω Ω∩Ω Ω∩Ω Ω∩Ω + ⋅ + ⋅ − ⋅ ⋅ + ∇ = ⋅ p p u p p p p p p p p p p p p p u p u p p p p d dS dS d V m d d V τ v x τ v u  u v v x σ v x b v x p δ δ α δ δ χ δ χ δ χ τ ( ) : & . (37) 42 where ∫ Ω∩Ω = p V d p p χ (x) x is the particle volume. Introducing a background grid and the grid shape function (x) i S that satisfies partition of unity 1 ) ( = Σi i S x , the admissible velocity field can be represented by the grid nodal data as =Σ i i i δv δv S (x) . Without the loss of generality, take u in Eq. (37) to be the displacement of the boundary particles at the current time step and u p u τ = σ n where u n is the unit outward normal to u ∂Ω . The momentum conservation, Eq. (37), can eventually be written for each node i as u i i b i i i p& = f int + f + f τ + f , (38) where the time rate change of nodal momentum =Σ / Δ , p i ip p p& S p t the nodal internal force vector int = −Σ ⋅∇ , p i p ip p f σ S V the nodal body force vector =Σ p p ip b i f m bS and the nodal traction force vector Σ ∫ ∂Ω ∩Ω = p i i p S dS τ f τ τ (x) . u i f is the force vector induced by the essential boundary condition given by Σ ∫ Σ ∫ ∂Ω ∩Ω ∂Ω ∩Ω = − p p i p p u i u i u p u p f σ n S (x)dS α (u  u)S (x)dS . (39) Sip is weighting function between particle p and node i given as ∫ Ω∩Ω = p S d V S p i p ip (x) (x) x 1 χ . (310) The weighting function in GIMP is C1 continuous and satisfies partition of unity. The momentum conservation, Eq. (38), can be solved at each node to update the nodal momentum, acceleration, and velocity. These updated nodal quantities can be 43 interpolated to the material particles to update the particles, as given by Bardenhagen and Kober (2004). It may be noted that the mass of each material particle does not change, so that the mass conservation equation is satisfied automatically. In the discretization of the weak form of the momentum conservation equation, a background grid is used. However, the computation is independent of the grid from one increment to another. Hence a spatially fixed structured grid can be used for convenience. In the background grid, no nodal connectivity is required and the integration is never performed on the element domain. Similar characteristics have been reported for other meshless methods, such as the meshless local PetrovGalerkin (MLPG) method (Atluri and Shen (2002)). For a uniform structured grid, the grid shape function in 3D is defined as the product of three nodal tent functions (Bardenhagen and Kober (2004)) S ( ) S (x)S ( y)S z (z) i y i x i i x = , (311) in which the nodal tent functions are in the same form, e.g., ⎪ ⎪ ⎩ ⎪ ⎪ ⎨ ⎧ ≤ − − − ≤ − ≤ + − − ≤ − ≤ − ≤ − = x i i x i x i x x i i x x i L x x x x L x x L x x L L x x x x L S x 0 1 ( ) / 0 1 ( ) / 0 0 ( ) . (312) Fig. 31 shows one 2D grid cell with four nodes. In this chapter, the particle characteristic function of the material particle located at (xp, yp) is taken as ( ) (x) y ( y) p x p p χ x = χ χ , (313) where ( x ) H[ x ( x l )] H[ x ( x l )] p x p x xp χ = − − − − + and H denotes the step function. 44 Fig. 31: 2D representation of a particle and a grid cell 3.3 Structured Mesh Refinement In this section, a refinement scheme for a structured mesh in GIMP is described. Since GIMP shares some characteristics with meshless methods, the GIMP method is expected to have h convergence in simulation (Atluri and Shen (2002)). The momentum conservation equation is essentially solved at each node (see Eq. (38)). Therefore, the number of equations to be solved is the same as the number of nodes. Finer grid and smaller material particles will lead to more accurate results. In some simulations, high stress gradients exist in small regions. For instance, in indentation with a sharp tip, the stress gradient is high in the workpiece beneath the indenter tip. In simulation of fracture problems, the stress gradient at the crack tip is high and of particular interest. Consequently, finer grid is needed for these regions; but away from these regions, a coarser grid can be used to reduce the computational cost. In conclusion, a uniform grid can be either too computationally expensive if it is too fine, or inaccurate, if it is coarse. A nonuniform grid with refinement can provide accurate results while minimizing the overall computational time. Grid refinement should maintain the same characteristics of the structured grid as much as possible, in order to replicate the grid generation in parallel processing. The proposed Lx Ly (x1 ,y1) (x2 ,y2) (xp ,yp) 2lx 2ly (x4 ,y4) (x3 ,y3) 45 refinement scheme is illustrated in Fig. 32 with one particle per cell assigned. The material particles that fill each cell are square in nature but denoted as circles for clarity. To understand this grid, one can consider that there are two overlapped structured grids. The coarse grid covers a rectangular region from (0, 0) to (8, 6) and the fine grid covers a region from (2, 2) to (6, 6). For each grid, the shape function can be evaluated from Eq. (312). To be consistent with any other general refinement, it is required that the region of the fine grid to be smaller than the coarse grid. Fig. 32: Refinement of structured grid with a refinement ratio of two When these two grids are merged into one, the shape function and the weighting function for the nodes at the boundary of the finer grid, for example, the nodes at (2, 2) and (2, 3) should be changed. These nodes are called transition nodes. To facilitate the computation of the interpolation function, define an influence zone for each node, denoted as [ − , + , − , + ] x x y y L L L L in 2D or [ − , + , − , + , − , + ] x x y y z z L L L L L L in 3D. The symbols in the square bracket define the size of the influence zone, whereas the subscript denotes the coordinate axis and the superscript denotes the direction of the axis. For example, − x L and + x L represent the sizes in the negative and positive X direction, respectively. The influence zone for each node in 2D is rectangular and it extends to the next immediate grid line to the left, X Y 0 2 4 6 8 4 2 6 46 right, bottom and top of the node. If no more grid lines exist in any direction, such as the boundary nodes, the size is zero in that direction. For example, in the refined grid in Fig. 32, the influence zones for the nodes at (2, 3) and (2, 4) are both [2, 1, 1, 1]. The influence zone for the node at (2, 2) is [2, 1, 2, 1]. Based on this definition, the influence zone for this node in the coarse grid is [2, 2, 2, 2], and in the uniformly fine grid is [0, 1, 0, 1]. Based on the influence zone, the nodal tent function in each direction can be modified as, for example, in the X direction, ⎪ ⎪ ⎩ ⎪ ⎪ ⎨ ⎧ ≤ − − − ≤ − ≤ + − − ≤ − ≤ − ≤ − = + + + − − − x i i x i x i x x i i x x i L x x x x L x x L x x L L x x x x L S x 0 1 ( ) / 0 1 ( ) / 0 0 ( ) . (314) Eq. (314) can be substituted into the grid shape function (Eq. (311)), and the weighting function between the particle p and the node i can be evaluated as ⎪ ⎪ ⎪ ⎪ ⎩ ⎪ ⎪ ⎪ ⎪ ⎨ ⎧ − − − ≥ − − − ≤ − + − ≤ − ≥ = − + + − − + otherwise l b a a L b L a l b a b a L b l b a b a L B L or A L S p x x p x p x x x x ip 2 /(2 ) /(2 ) 0 2 ( ) /(2 ) 0 2 ( ) /(2 ) 0 2 2 2 2 2 2 . (315) where i p A = x − x −l , i p B = x − x + l , a =max( , − − ) x A L and b =min( , + ) x B L . When − = +x x L L , Eqs. (314) and (15) are degraded to the cases for uniform grid. Without detailed proof, the grid shape function and the weighting function still satisfy partition of unity. Similarly, the gradient of the modified weighting function can be computed. 47 It may be noted that in Fig. 32 the refinement ratio is two, i.e., the length of a side of a coarse cell is twice that of the fine cell. To maintain the convenience of the structured grid, only integer refinement ratio should be used. All nodal positions can be computed from the domain of each grid and the cell sizes. The proposed refinement scheme can be applied to any integer refinement ratio and for multiple times for successive refinements. As an example, the weighting function between a particle of size 0.5×0.5 and a node at (0, 0) with an influence zone of [1, 1, 1, 1] is shown in Fig. 33 (a). The particle is on the XY plane and the weighting function is computed at each particle position. For comparison, the influence zone is changed to [1, 0.5, 1, 1], representing a transitional node, while other conditions are the same. The weighting function for this case is plotted in Fig. 33 (b). It can be seen that the weighting function for the transitional node is still C1 continuous. Fig. 33: Effect of influence zone on the weighting functions 3.4 Numerical Simulations 3.4.1 Tracking particle deformation Prior to presenting the results of structured mesh refinement, tracking particle deformation is addressed first since this is necessary in later simulations to achieve accuracy. The material particles are initialized into regular shapes, normally square and (a) Influence zone: [1, 1, 1, 1] (b) Influence zone: [1, 0.5, 1, 1] 1.5 1 0.5 0 0.5 1 1.5 1 0.5 0 0.5 1 1.5 0 0.2 0.4 0.6 0.8 Y Weight X 1.5 1 0.5 0 0.5 1 1.5 1 0.5 0 0.5 1 1.5 0 0.2 0.4 0.6 0.8 Y Weight X 48 cube for 2D and 3D simulations, respectively. All the physical quantities in a particle domain are considered to be uniform. The shape of the particle changes during deformation. So, it is important to track the deformation of each particle. Fig. 34 illustrates the deformation of the particles in 2D when the particles are stretched in the Xdirection. If the particle deformation is not tracked, gaps will form between neighboring particles in the Xdirection, as shown in Fig. 34 (a). Due to the Poisson’s effect, there will be overlapping between particles in the Ydirection, if the particles do not follow the deformation of the materials properly. When the stretch and gaps are large enough, the particles would be separated. Fig. 34 (b) shows the correct deformation in which contiguous particles remain contiguous after deformation. Fig. 34: Schematics showing overlaps and gaps that may occur when particle deformation is not tracked To track the particle deformation, a convenient approach is to calculate the deformed particle shape based on strain history. Since the linear strain increment is computed in GIMP, the effectiveness and validity of tracking particle deformation based on strain would be limited to relatively small deformations. As an example, Fig. 35 (a) shows the simulation of a uniaxial tension problem under plane strain conditions. The material is silicon. Its mass density is 2.71 g/cm3, Young’s modulus 175.8 GPa, and the Poisson’s Overlap X Y (a) No particle deformation (b) Actual particle deformation Gap 49 ratio 0.28. The background grid size is 0.002 × 0.002 mm2. One particle per cell is assigned initially and the time step is 0.02 ns (1 ns = 109 s). The applied pressure increases linearly with time from 0 to 10 ns and is then maintained constant, as shown in Fig. 35 (b). The elongation in the Xdirection of the particle is computed as (1 ) 0x x +ε l , where 0 x l is the initial length of the particles. Separation of the particles occurred at ~25% strain at 6 ns before the full pressure was applied, as shown in Fig. 35 (c). Similar problems have been reported by Guilkey (2005). Fig. 35: Simulation showing separation when the particle deformation is tracked by strain Tracking particle deformation by strain could be more effective if nonlinear strain is used in the GIMP method. However, recovery of the deformed shape based on nonlinear strain involves additional complications. Based on the GIMP algorithm, there is another convenient approach to track the particle deformation. Numerically, the displacement and velocity of each particle in GIMP are computed at the center of the particle to represent X Y L=0.06 mm H=0.04 (a) (c) p (b) t p 40 GPa 0 10 ns X Y 0 0.02 0.04 0.06 0.08 0.1 0.01 0.02 0.03 0.04 0.05 σx 95000 85500 76000 66500 57000 47500 38000 28500 19000 9500 0 Separation (MPa) 50 the entire particle domain. However, in reality, the velocity and deformation at the corners of a particle can be different from the center. Fig. 36 shows four 2D contiguous particles sharing one common corner point at the middle. This corner point should have unique displacement and velocity. As a result, it is helpful to track the displacement and velocity of each corner to track the particle deformation. It is not difficult to compute the velocity of the particle corner given its location. It is computed from the interpolation from the background grid, similar to the center of the particle. For a 2D particle, in addition to updating the position of the center of the particle, the positions of the four corners are updated at each increment. To compute the weighting function for the corners, a fictitious size can be assigned to each corner. Numerical simulation shows that the result is not sensitive to this size in the range of 10% to 80% of the initial particle size. The new particle shape can be obtained by connecting the four corners with straight lines. In order to avoid numerical integration of the interpolation function, it is assumed that the deformed material particle shape is rectangular with edges parallel to the coordinate axes. The size of the rectangle is, therefore, determined from the extent of the corners. As will be demonstrated later in this section, this assumption does not introduce any significant error while it can greatly improve the efficiency of the GIMP algorithm. Fig. 36: Velocity field of a continuum region (the arrows represent both direction and magnitude) 1 2 4 3 51 Fig. 37: GIMP results with tracking deformation of corners and their comparison with FEM Using this approach to track the particle deformation, the problem in Fig. 35 (a) was simulated again with GIMP and the results at 20 ns are plotted in Fig. 37 (a). No separation of particles was seen during the entire simulation up to 50% strain. It is noted that each material particle is plotted as a square of the same size and the particle deformation is not shown due to software limitations on visualization. For comparison, the same problem is simulated using FEM (Abaqus/Explicit) and the FE result is shown in Fig. 37 (b). It can be seen that these two sets of results agree reasonably well with each other; the maximum difference in maximum tensile stress is ~8%. 3.4.2 Indentation problem To verify the refinement algorithm, a 2D indentation problem was simulated. Pressure is applied at the top of the workpiece, as shown in Fig. 38. The dashed lines indicate the (b) FEM (MPa) (a) GIMP 52 borders of refinement levels. The work material is silicon with the properties the same as those given in the previous section. Due to symmetry with respect to the Yaxis, only half of the workpiece is modeled. The size of the model is 0.027 × 0.04 mm2. The pressure increases linearly with time from t = 0  30 ns and is then kept constant at 60 GPa. Fig. 38: Twodimensional indentation simulation Several simulations were performed under different settings for the purpose of comparison. In the first simulation, a uniform grid with a cell size of 0.001 × 0.001 mm2 is used and the time step is 20 ps (1 ps = 10–12 s). The stress distribution in the workpiece at t = 20 ns is shown in Fig. 39 (a). In this figure, the units of length and stress are mm and MPa, respectively. In the second simulation, as indicated in Fig. 39 (b), two levels of refinements are used with the refinement ratio to be 2. The cell lengths are 0.002 mm and 0.001 mm for the first and second levels, respectively. The fine level covers a rectangular area of the workpiece from (0, 0.02) to (0.02, 0.04). The grid is fixed in space; the material particles initially in the fine region move to the coarse region during deformation. As shown in Fig. 39 (b), some fine material particles have moved below the line Y = 0.02 mm. It may be noted that in Fig. 39 the material particles are plotted as squares corresponding to their initial sizes. Gaps between particles are intentionally shown to depict the particle sizes. X Y p Level 2 Level 3 Symmetry p t 60 GPa 0 30 ns Level 1 53 Fig. 39: Comparison of the stress distributions at different levels of refinements in force indentation Three levels of refinements are used in the third simulation of the same problem. In this simulation, the time step is 10 ps and the results are shown in Fig. 39 (c). The stresses agree very well with the previous two simulations. Fig. 39 (d) shows the result when the particle deformation is not tracked. A severe material separation was observed at t = 36 ns, as indicated by the arrow. Additionally, the displacement history of the particle in the middle of the top surface for the third simulations is shown in Fig. 310. The same result of the integration point of the element in the middle of the top surface in FE simulation using ABAQUS/Explicit is also shown in Fig. 310 for comparison. It can be seen that Separation (a) One level (b) Two levels (c) Three levels (d) Deformation not tracked (MPa) 54 the displacement compares well with FE up to t = ~20 ns. After this time, FE simulation aborted due to mesh distortion. This demonstrates the capability of GIMP using a grid with structured refinement in handling large deformations. Time (ns) Uy (μm) 0 10 20 30 40 8 6 4 2 0 GIMP FEM FEM aborted Fig. 310: Comparison of displacement history with FE It may be noted that if the exterior corners of the surface particles are tracked from the nodal interpolations in the same way as the interior corners, simulations tend to become unstable due to erroneous surface corner displacements. This problem was observed to be strictly and consistently associated with the particles with external tractions applied. It is caused by insufficient nodal interpolation. To eliminate this problem, the exterior corners of the surface particles were tracked by strain only, as used in these simulations. 3.4.3 Verification of the displacement boundary condition Next, the results on the validation of the displacement boundary conditions are presented. For dynamic simulations, an artificial damping may be introduced. With damping, the nodal momentum can be updated as d t i u i i b i i i Δp = (f int + f + f τ + f − p )Δ , (316) where d is the artificial damping coefficient. 55 A rectangular slab is fixed on the left and a displacement boundary condition is applied on the right, as shown in Fig. 311. The material is silicon and its properties are given in the previous section. The cell size in this simulation is 0.5 × 0.5 mm2 and four particles are assigned to each cell initially. The time step is 5 ns. The prescribed displacement increases linearly with time to 0.5 mm at t = 10 μs, which corresponds to a velocity of 50 m/s, and remains constant thereafter. This problem is simulated in FE using Abaqus/Explicit for comparison. The displacement in the Xdirection, Ux, for a particle initially centered at (2.625, 1.625) as a function of time is plotted in Fig. 312 (a). It can be seen that without damping, the vibrations of displacement from both FEM and GIMP simulations are in phase before t = 14 μs, but out of phase afterwards. The dashed line represents the steady state displacement at this point. It can be seen that when an artificial damping of 106 s1 is used, the GIMP solution converges quickly. The error of the converged displacement is 0.46%. Fig. 312 (b) shows the comparison of the normal stress in the Xdirection. Good comparison between GIMP and FE has been obtained. With damping, the GIMP solution converges to the analytical solution for static simulations. Fig. 311: A simple tension problem with displacement boundary conditions X Y L=6 mm H=3 mm u 56 Fig. 312: Comparisons of the displacement and stress 3.4.4 Stress concentration problem Fig. 313 shows a copper plate (60 × 60 mm2) with a central hole (4 mm diameter) subjected to a distributed load. This problem is simulated using GIMP as a dynamic problem with a damping factor of 1000 s1 and three levels of refinement. The cell sizes at these three levels are 1.0 mm, 0.5 mm, and 0.25 mm, respectively. One particle is assigned to each cell not adjacent to the hole initially. The cells close to the circular hole are assigned 25 particles each to model the circular edge more accurately with the use of square particles. It may be noted that all the particles occupy square areas initially. The time step is 10 ns and the applied distributed load intensity p = 10 MPa. Fig. 313: A plate with a circular hole subjected to tension The stress distribution after 4000 steps is shown in Fig. 314 (a) and the area close to the hole (bottom left corner in Fig. 314 (a)) is magnified in Fig. 314 (b). The normal stress p p X Y θ Time (μs) Ux (mm) 0 10 20 30 40 0 0.05 0.1 0.15 0.2 0.25 0.3 GIMP d=0 FE d=0 GIMP d=106 1/s Steady State (a) Displacement history Time (μs) σx 0 10 20 30 40 0 5 10 15 20 GIMP d=0 FE d=0 GIMP d=106 1/s Steady State (GPa) (b) Stress history 57 of the particle at the top of the hole is 30.7 MPa when the applied tension is 10 MPa. This gives a stress concentration factor of 3.07. The tangential stress of the particles along the hole circumference, normalized by the applied pressure, is plotted in Fig. 315 in comparison with the theoretical value (Pilkey (1997)). A good agreement between the numerical simulation and theoretical value is obtained. This demonstrates that the GIMP refinement algorithm is effective for problems involving significant stress gradients. Furthermore, with the use of small square particles, GIMP is capable of modeling curved surfaces . Fig. 314: Normal stress in the Xdirection with p = 10 MPa Angle θ (degrees) Normalized Tangential Stress 20 15 30 45 60 75 90 1 0 1 2 3 4 Numerical Theoretical (Pilkey, 1997) Fig. 315: Normalized tangential stress along the circumference of the hole (a) Overall (b) Magnified (MPa) 58 3.4.5 Static stress intensity factor Next, the GIMP refinement algorithm is implemented to determine the stress field and stress intensity factor for a mode I crack problem to determine the capability of the GIMP refinement algorithm in simulation of stress distribution near a crack. Guo and Nairn (2003) have successfully extended the MPM method to compute stress distribution in a plate with explicit cracks using multiple nodal fields along the crack surface. The physical quantities of material particles on each side of the crack were interpolated using variables in the field on that side of crack surface. In their simulation, a uniform mesh was used. Since the stress gradient at the crack tip is very high, a refined mesh near the crack tip and a coarse mesh in the far field should lead to savings in computational time while maintaining the same accuracy with the use of a uniform fine mesh. The fracture problem is thus an appropriate problem to evaluate the refined GIMP algorithm. For this purpose, the same fracture problem by Guo and Nairn (2003) using MPM is modeled, that is, a double cantilever beam (DCB) with a crack as shown in Fig. 316. Fig. 316: Geometry of a double cantilever beam with a crack In the area close to the crack tip, finer meshes are used, while in the area far away from the crack tip, coarse meshes are used. The thickness of the plate is 1 mm, thereby justifying a plane stress condition for this problem. The material of the DCB is L=10 2H=2 a=5 X Y F F Unit: mm Crack 59 considered to be homogeneous, isotropic, and linearly elastic. Its density, Young’s modulus, and Poisson’s ratio are 1500 kg/m3, 2300 GPa, and 0.33, respectively. The applied force is F = 4×10−4 N and this results in a mode I crack problem. The static stress intensity factor for the DCB can be calculated using the following equation (Kanninen 1973) 3 / 2 2 3 ( 3 / 2) H K F a H I + = . (317) This problem has been simulated using MPM with uniform grids of three sizes, 4 mm, 2 mm, and 1 mm (Guo and Nairn (2003)). They have computed the energy release rate and determined the stress intensity factor from Jintegral. Their results indicate that the stress intensity factor determined from finer grid is more accurate and closer to the theoretical value. In the mesh refinement GIMP algorithm used in this study, the energy release rate was computed using the virtual crack closure technique (Rybicki and Kanninen (1977); Wang, Karuppiah, Lu, Roy, and Komanduri (2005)). The energy released during an infinitesimal crack growth of Δa is assumed to be the energy required to close the crack to its initial size. Hence, the energy release rate G is determined by ∫Δ Δ → ⋅ Δ Δ = a a da a G 0 0 ( ) 2 lim 1 σn u , (318) where n is the crack surface normal. For the 2D mode I crack in the Xdirection, ∫Δ Δ → ⋅Δ Δ = a I a yy yσ u da a G 0 2 0 lim 1 . (319) In GIMP, the energy release rate can be computed as 60 ( ) 2 1 1 2 I tip y y F u u t a G − Δ = , (320) Fig. 317: Material points and background grid around the crack tip where the superscripts 1 and 2 denote the two material particles immediately to the left of the crack tip as shown in Fig. 317, Δa is the X distance between particle 1 and the crack tip, t is the thickness of the beam. tip F is the nodal force to hold the crack tip together (Rybicki and Kanninen (1977)) and is computed as the crack tip nodal force from one side of the crack in this simulation. The mode I stress intensity factor for the static crack is given by K GE I = . (321) GIMP simulations were carried out using four material points per cell. The time step is 0.1 μs, and a damping coefficient of 4000 s1 is used to allow the results to converge to the static data. The computed stress intensity factors using uniform grids of 4 mm and 1 mm, respectively, are plotted in Fig. 318 (a) and the theoretical value calculated from Eq. (317) is also shown for comparison. For simulation with refinement, two levels of refinement were used, i.e., 2 mm grid size for the coarse level and 1 mm grid size for the fine level, and the extent of the fine level is indicated by the dashed square in Fig. 316. The computed stress intensity factor using the refinement algorithm is also plotted in Fig. 318 (a). It is seen that the results from two levels of refinement are identical to the Crack surface X Y 1 2 61 results from one uniform fine level with 1 mm cell size. Moreover, the simulation time using the refinement algorithm is 38% shorter than that of one uniform fine level. The computed stress intensity factor became even closer to the theoretical value when a third level of refinement was added for the crack tip as shown in Fig. 318 (b). For the case of using uniform grid of 2 mm, the computed stress intensity factor using GIMP and MPM (Guo and Nairn (2003)) are plotted in Fig. 319. The MPM results were computed using a damping factor of 1000 s1, and therefore, more oscillations can be seen as expected. Fig. 318: Computed stress intensity factor as a function of time Fig. 319: Comparison of the stress intensity factor with MPM results (a) Overall (b) Magnified 62 Fig. 320: Stress distribution in the beam using three levels of refinements at t = 2.5 ms Fig. 320 shows the distribution of σy at t = 2.5 ms when the force applied on the beam was changed to F = 4 N, 10000 times of the previous value, with other parameters the same. In this simulation, three levels of refinements with cell sizes of 1 mm, 0.5 mm and 0.25 mm, were used and the different density of material points in each level can be clearly seen in the figure. The computed stress intensity factor, scaled by 10000 times, still compares well with the theoretical value. It is noted that deformation near the crack surfaces has caused material points crossing cell boundaries, a situation where MPM would give numerical noise such as alternating stress signs. Despite the extent of the deformation, and material points crossing cell boundaries, GIMP method with the use of particle deformation tracking still gives correct results, further verifying the structured refinement methods developed herein. 3.5 Conclusions 1. A spatial refinement scheme for a structured grid was developed by adding transitional nodes and by changing the influence zone of the transition nodes in GIMP. The influence zone is square for uniform grid nodes and rectangular for transitional nodes. This influence zone affects the computation of the nodal shape function. The computation of the weighting function remains the same as for the uniform grid. The refinement scheme (MPa) Crack 63 can be applied successively and the refined grid remains structured in each refinement level, i.e., every node can be determined by the extent of the grid level and cell size. The refinement scheme was implemented and several problems such as tension, indentation, stress concentration, and stress distribution near a crack (mode I crack problem) were modeled to demonstrate its effectiveness and accuracy. A good agreement has been obtained between numerical and theoretical results, indicating the validity of the structured mesh refinement for GIMP scheme. 2. The GIMP algorithm has also been extended to include the displacement boundary condition, based on the approach used in the meshless local PetrovGalerkin (MLPG) method, Atluri and Zhu (2000). A penalty parameter is used to impose the displacement boundary condition and a nodal force vector because the displacement boundary condition is introduced to the nodal momentum governing equation. A uniaxial tension problem with constant pulling velocity was simulated to verify the displacement boundary condition. 3. A method to track the material particle deformation was developed and verified in one example. When the particle deformation is not tracked, artificial separation was observed when the particle strain increases to a certain level. In tensile simulations, when the normal strain is ~25%, material particles tend to separate from the body. Our approach tracks the displacement of each corner of the material particles. Since neighboring particles share corners, no separation would occur during deformation using this approach. 64 Chapter 4 Simulation of Deformation Mechanism and Failure of Bulk Metallic Glass (BMG) Foam Using the GIMP Method A new bulk metallic glass (BMG) foam, processed by thermoplastic expansion, has recently been discovered (Schroers et al. (2003); Brothers and Dunand (2004); Hanan et al. (2005)). In this investigation, the microstructure and mechanical properties of BMG foam were determined using microtomography combined with insitu compression. Numerical simulation of the mechanisms of deformation and failure of the BMG foam, which has never been done before, is conducted in this investigation. The complex cellular geometry of the foam was modeled using the generalized interpolation material point (GIMP) method by creating material points based on grayscales at voxels determined from microtomography. Nanoindentation tests were conducted to determine the local elastic modulus of the foam matrix and the results are used as the input to the GIMP simulation. Several cubic volume elements (from 0.5763 mm3 to 1.7283 mm3) taken from the 3D tomographic image were used in the simulation. The size of the cube is increased until further increase in volume does not lead to a change in mechanical response; the cube with this size is determined as the representative volume element. The simulated compressive stressstrain curve using GIMP was compared with the experimental data and a good agreement on the initial slope (representative of the average Young’s modulus of the foam) was obtained. The stress contour in the foam was 65 illustrated by analyzing different sections of the foam at different simulation times. The densification and hardening of the ductile metallic foam under compression was simulated to demonstrate the capabilities of GIMP and to study the failure mechanisms. 4.1 Introduction Metallic glass foam has been used in various engineering applications due to such combination of unique properties as low density, high specific strength, high energy absorption, and superior thermal insulation. Foam materials are generally divided into opencell and closedcell foams. In general, opencell foams have higher porosity and softer behavior than closedcell foams. The uniaxial mechanical behavior of a foam is characterized by three regimes, namely, initial elastic response, the stress plateau that follows, and final stiffening due to densification (Gibson and Ashby (1997); Miller (2000); Katti et al. (2006)). Theoretical modeling of the foam mechanics has been carried out for both regular and irregular opencell foams. The nonlinear stressstrain relation of the foams at high strains was modeled by the buckling of elastic cell edges (Zhu et al. (1997)) or by the implementation of nonlinear kinematics (Wang and Cuitino (2000)). Most numerical simulations of foam materials reported in the literature were either at the macroscopic scale or the mesomechanical (micromechanical) scale. At the macroscopic scale, the foam material is considered as a homogeneous material with homogenized properties determined from measurements such as uniaxial compressions (Gibson (1997)). The work using this approach is rather limited, due primarily to the wellknown problem associated with simulating cells in contact. In this approach, large foam structures can be simulated at relatively low computational cost and the deformations under different loading conditions can be predicted reasonably well (Meguid et al. (2002); Meguid et al. 66 (2004)). However, such simulations require knowledge of loading conditions a priori and the associated material model may not be readily available in some situations (Wicklein et al. (2004)). In the mesomechanical simulations, the microstructure of the foam is modeled explicitly. Xray tomography is used to determine the threedimensional (3D) microstructure and the 3D geometry is used in a numerical model for simulation. In this approach, the discrete element/particle size is dictated by the resolution of the Xray tomography. Kadar et al. (2004) investigated the mechanical behavior of closedcell Al porous foam by indentation using the finite element method (FEM). Wicklein and Thoma (2005) conducted a series of FEM simulations and determined the elastic and plastic response of an opencell aluminum foam. The mesomechanical approach provides microstructural evolution during deformation; however, the computational cost can be very high. For model generation using mesomechanical simulation approach, Wicklein et al. (2004) described three possible approaches to discretize a foam structure into FEM meshes. In the first approach, beam and plate elements are constructed based on the cell votices in the tomographic images. In the second approach, voxels are transformed into cubic elements in a predefined structured mesh. In the third approach, tetrahedral elements are used to mesh the entire foam structure. The second approach was used successfully for the simulation of aluminum foams (Wicklein et al. (2004); Wicklein and Thoma (2005)) under a few percent deformations. However, simulation of a fully densified foam is a challenge in FEM because of large distortions of the microstructure involved and the internal contact of cell walls upon closure of cells. Brown et al. (2000) used hexahedral elements to model the struts of an opencell foam under impact loading. It was pointed 67 out that the drawbacks of this approach involve the use of a large number of elements and small timesteps. It is noted that the struts all have regular cross sections so that it was relatively easy to mesh all the struts. Recently, the generalized interpolation material point (GIMP) method was introduced for dynamic simulation of materials (Bardenhagen and Kober (2004); Ma et al. (2005)). GIMP has evolved from the material point method (MPM), originally developed by Sulsky, Zhou and Schreyer (1995) and subsequently refined by others (e.g., Tan and Nairn (2002); Wang et al. (2005)). In GIMP, the workpiece is discretized into a collection of material particles (points). The material points carry all the physical variables, including the mass, momentum and stress, in the simulation. A background grid, usually a structured grid fixed in space, is used to discretize the momentum conservation equation and solve it at the grid nodes using explicit time integration. In a regular GIMP computational step, the material point information is first interpolated to the grid nodes. Then, the momentum conservation equation is solved at the nodes to update the nodal information. Next, the updated nodal information is interpolated to update the position, velocity, strain, and stress at the material points. Finally, the background grid is redefined (usually reset the nodal quantities to zero but keep the nodal positions fixed). GIMP shares some similarities with other meshless methods, e.g., Atluri and Zhu (2000). Fig. 41 (a) illustrates an elliptical workpiece discretized into a collection of material points of different sizes with the background grid shown. When two bodies are close enough, as indicated by open and closed circles in Fig. 41 (b), the material point information will be interpolated to the same grid node in the middle when one background is used. In this case, these two workpieces are in nonslip contact. The 68 intrinsic features associated with GIMP, such as ease of discretization, no bodyfixed meshes used so that mesh distortion is not an issue and natural capability of handling nonslip contact, made the simulation of some complex problems possible. Fig. 41: Illustrations of the discretization scheme and intrinsic contact in GIMP GIMP was used successfully to simulate densification of opencelled foam materials (Bardenhagen et al. (2005); Brydon et al. (2005)). In GIMP simulations, each voxel in the Xray tomograph is converted into a material particle and an Eulerian structured grid is used for solving the field equations. A representative volume element (RVE) was used to simulate the bulk response of a foam. Their results have indicated that the apparent Poisson’s ratio is negative in the stress plateau and the dynamic response of the compression is dominated by a compression wave, which is an inertial effect absent in homogeneous material and its velocity is much slower than any characteristic wave speeds (Brydon et al. (2005)). In this work, we report the results of a simulation of a closedcell foam made of bulk metallic glass. Most metal and its alloys in engineering applications are in the form of crystalline structures; their mechanical properties are dependent on the crystal structure, chemistry, microstructure, and processing conditions. Metallic glass is another form of a metal or an alloy. In a metallic glass, the material is amorphous and has no longrange atomic order. Fig. 42 (a) and (b) are TEM images of the microstructure of a crystalline (a) Discretization and background grid (b) Nonslip contact 69 zirconium alloy and the same alloy in the amorphous form (Hufnagel (2005)). Metallic glass can be formed by rapid cooling of a liquid metal; the cooling rate is so high that ordered crystalline structure cannot be developed in the metal, resulting in amorphous structures. Fig. 42: High resolution TEM images of (a) crystalline and (b) amorphous structures of a zirconium alloy (Hufnagel (2005)) Metallic glass foams recently have attracted great attention due to their unique properties. They exhibit high specific strength, high stiffness, high energy absorption (Ashby et al. (2000)), and more importantly they can be used under those conditions, including, high temperatures where other foams (e.g., polymer foams) cannot be used. BMG foams combine exceptional strength, elasticity, wear, and corrosion resistance with modest densities and low processing temperatures. They are being considered for structural and biocompatible implant applications (Hanan et al. (2005); Brother and Dunand (2005)). The thickness of BMG material is usually limited (on the order of a few millimeters to a centimeter or so) by the requirement that a high enough cooling rate in the center of the wall is necessary to form the amorphous structures. For a 



A 

B 

C 

D 

E 

F 

I 

J 

K 

L 

O 

P 

R 

S 

T 

U 

V 

W 


