MULTISCALE MODELING OF POLYMER NANOCOMPOSITES By Azadeh Sheidaei A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of Mechanical Engineering - Doctor of Philosophy 2015 ABSTRACT MULTISCALE MODELING OF POLYMER NANOCOMPOSITES By Azadeh Sheidaei In recent years, polymer nano -composites (PNCs) have increasingly gained more attention due to their improved mechanical, barrier, thermal, optical, electrical and biodegradable properties in comparison with the conventional micro -composites or pristine polymer. With a modest addition of nanoparticles (usually less than 5wt. %), PNCs offer a wide range of improvem ents in moduli, strength, heat resistance, biodegradability, as well as decrease in gas permeability and flammability . Although PNCs offer enormous opportunities to design novel material systems, development of an effective numerical modeling approach to p redict their properties based on their complex multi -phase and multiscale structure is still at an early stage. Developing a computational framework to predict the mechanical properties of PNC is the focus of this dissertation. A computational framework ha s been developed to predict mechanical properties of polymer nano -composites. In chapter 1, A microstructure inspired material model has been developed based on statistical technique and this technique has been used to reconstruct the microstructure of Hal loysite nanotube (HNT) polypropylene composite and exfoliated Graphene nanoplatelet (xGnP) polymer composite . Chapter 2 is the summary of the experimental work to support the numerical work. Melt extrusion followed by injection molding was used to manufacture high density polyethylene (HDPE) Œ xGnP nanocomposties. Scanning electron microscopy (SEM) also was performed to determine particle size and distribution and to examine fracture surfaces. Particle size was measured from these images and has been used for calculating the probability density function for GNPs in chapter 1. A series of n anoindentation tests have been conducted to reveal the spatial variation of the superstructure developed al ong and across the flow direction of injection -molded HDPE/GNP. The uniaxial tensile test and shear test have been conducted on HDPE and xGnP/HDPE specimens. The stress -strain curves for HDPE obtained from these experiments have been used in chapter 5 to c alibrate the modified Gurson ŒTvergaard ŒNeedleman to capture the damage progression in HDPE. In chapter 3, the 3D microstructure model developed in chapter 1 was incorporated in a damage modeling problem in nanocomposite where damage initiation has been mo deled using cohesive -zone model . The finite element model of progressive debonding in nano -reinforced composite has been proposed b ased on the cohesive -zone model of the interface. In order to model cohesive -zone , a cohesive zone traction displacement relation is needed. Thi s curve may be obtained either through a fiber pullout experiment or by simulating the test using molecular dynamics. In the case of nano -fillers, conducting fiber pullout test is very difficul t and result is often not re producible . In chapter 4, molecular dynamics simulation of polymer nanocomposite has been performed. One of the goals was to extract the load -displacement curves of graphene/HDPE pullout test and obtain cohesive zone parameters in chapter 3. Finally, in chapter 5, a damage model of HDPE/GNP nanocomposite has been developed based on matrix cracking and fiber debonding . This 3D microstructure model was incorporated in a damage modeling problem in nanocomposite where damage initiation and damage progression have been modeled using cohesive -zone and modified Gurson -Tvergaard -Needleman (GTN) material models. iv This dissertation is dedicated to my parents (Soraya Sheidaei & Nasser Sheidaei ) for their love, endless support and encouragement. v ACKNOWLEDGEMENT S I would like to express my gratitude to all those who gave me the possibility to complete this dissertation . I want to express my sincere gratitude to my advisor, Professor Farhang Pourboghrat , who throughout my doctoral studies has contributed with excellent support and encouragement to commence and achieve this work. I have been amazingly fortuna te to have an advisor who gave me the freedom to explore on my own . I have to thank Prof. Lawrence Drzal who supported me with excellent scientific help, particularly in the domain of processing polymer nanocomposites . I am grateful to the CMSC staff at M SU, Mike Rich, Brian Rook, Dr. Per Aske land and Ana Veskovic . Their generous help and assistance enable d me to set up and run experiments smoot hly. I also want to express my gratitude to Prof. Thomas Pence and Prof. Xinran Xiao for serving on my PhD commi ttee, for their interest in my research and careful evaluation of this dissertation. I wish also to express my deepest gratitude to Dr. Majid Baniassadi from University of Strasbourg in France for his valuable ideas and suggestions and fruitful discussions. His encouragements have been a major reason for me to start and advance this work. Special thanks are reserved for Dr. Taejoon Park from Michigan State University and Dr. Fadi Abu -Farha from Clemson University. I am also grateful to the Compo site Vehicle Research Center for the financial support . Finally, I am heartily thankful to my dear mentor, friend, soulmate and husband, Dr. Shahram Pouya , who has always been there for me, understanding and unconditionally supportive of my endeavors as I pursued this goal. vi TABLE OF CONTENTS LIST OF TABLES ..................................................................................................................................... viii LIST OF FIGURES ..................................................................................................................................... ix KEY TO SYMBOLS .................................................................................................................................. xiii Chapter 1 3-D microstructure reconstruction of polymer nano -composite using FIB ŒSEM and statistical correlation function . Application to halloysite clay nanotubes polypropylene composite ............................................................................................... 1 1.1 Introduction and Literature Review ............................................................................... 1 1.2 Statistical correlation function ....................................................................................... 4 1.3 Materials and synthesis ................................................................................................. 7 1.4 Mechanical properties .................................................................................................... 7 1.5 Serial sectioning of the nano -composite using FIB ŒSEM ............................................... 8 1.6 3-D r econstruction using VCAT® software .................................................................. 11 1.7 3-D reconstruction of nano -composite using statistical correlation functions ............ 13 1.7.1 Reconstruction procedure ............................................................................................ 14 1.8 Finite element mesh using vcat2tets software ........................................................... 17 1.9 Finite element analysis of RVEs .................................................................................... 18 1.10 3D Reconstruction of xGnP/Polymer Using Statistical Correlation Functions ............. 21 1.11 Conclusion and future work ........................................................................................ 24 REFERENCES ................................................................................................................. 26 Chapter 2 Processing and testing of nanocomposites .................................................................. 31 2.1 Processing of nanocomposites ..................................................................................... 31 2.2 Processing of Polymer/Graphene Nanocomposites .................................................... 32 2.3 Processing of high density polyethylene -exfoliated graphene nanoplatelet nanocomposites .......................................................................................................... 34 2.4 Scanning Electron Microscopy ..................................................................................... 38 2.5 Nanoindentation .......................................................................................................... 40 2.6 Mechanical Properties .................................................................................................. 45 2.7 Shear Test ..................................................................................................................... 48 REFERENCES ................................................................................................................. 51 Chapter 3 An interfacial debonding -induced damage model for graphite nanoplatelet polymer composites ................................................................................................................... 54 3.1 Introduction and literature survey ............................................................................... 54 3.2 Nanocomposite model ................................................................................................. 58 3.2.1 Representative volume element (RVE) ........................................................................ 58 3.2.2 Cohesive zone model .................................................................................................... 59 3.3 Material behavior and FE simulation ........................................................................... 65 3.4 Results and discussion ................................................................................................. 67 3.4.1 The effect of the GNP™s volume fraction and aspect ratio in perfectly bonded nanocomposite ............................................................................................................. 67 vii 3.4.2 Comparing the effect of the GNP™s volume fraction and aspect r atio in perfectly bonded and cohesively bonded nanocomposites ........................................................ 69 3.4.3 The effect of the GNP™s aspect ratio and volume fraction in weakly bonded nanocomposite ............................................................................................................. 73 3.5 Summary and conclusions ............................................................................................ 77 REFERENCES ................................................................................................................. 79 Chapter 4 Molecular dynamics simulation of polymer nanocomposites ..................................... 83 4.1 Introduction .................................................................................................................. 83 4.2 Theory ........................................................................................................................... 84 4.3 United -atom model ...................................................................................................... 84 4.4 Fundamental algorithms in MD simulation .................................................................. 87 4.5 Ensemble average in MD simulation ............................................................................ 89 4.6 Nanoscale RVE .............................................................................................................. 90 4.7 Measuring interphase thickness for Graphene/PE system .......................................... 93 4.8 Comparing interfacial properties of CNT and Graphene in polymer ........................... 94 4.9 Conclusion .................................................................................................................... 97 REFERENCES ................................................................................................................. 99 Chapter 5 Damage model for thermoplastic polymer nanocomposite ...................................... 101 5.1 Damage model for high density polyethylene ........................................................... 102 5.1.1 Modified Gurson -Tvergaard -Needleman (GTN) model .............................................. 103 5.1.2 Calibration of modified GTN model for high density polyethylene ........................... 114 5.2 Damage model for HDPE/GNP nanocomposite ......................................................... 125 5.2.1 Representative volume element ................................................................................ 125 5.2.2 Boundary Conditions .................................................................................................. 126 5.2.3 Results ........................................................................................................................ 128 5.3 Conclusion and future work ....................................................................................... 129 REFERENCES ............................................................................................................... 132 viii LIST OF TABLES Table 1.1 . Properties of Halloysite clay nanotube and polypropylene (Lu, Chen et al. 2011) (top) and, elastic constants measured using real RVE and statistical RVE (bottom) .. . 21 Table 2.1 . Geometrical and surface characteristics of GNP -15 (Jiang and Drzal 2010) .. .............. 35 Table 2.2 . Elastic modulus measured in the middle of nanocomposite . ..................................... 44 Table 2.3 . Elastic modulus measured near the edge of nanocomposite . ..................................... 45 Table 3.1 . Cohesive zone model parameters for opening and sliding modes . ............................ 65 Table 3.2 . Cohesi ve parameters for weak bonding . ..................................................................... 74 Table 4.1 . Parameters for the bond stretching, the bond angle bending, and the torsional or dihedral potentials to model the intramolecular interaction of a linear polyethylene chain in the framework of an united -atom approach (Hossain, Tschopp et al. 201 0). .................................................................................................... 86 Table 5.1 . Properties of K46 -06-185 (High density polyethylene copolymer) . .......................... 102 Table 5.2 . Hardening evolution of the undamaged material . ..................................................... 122 Table 5.3 . Material coefficients for modified G -T-N model . ....................................................... 122 ix LIST OF FIGURES Figure 1.1 . Four different configurations of vectors for two -point correlation function . ............... 6 Figure 1.2 . Tensile modulus of HNT/PP and PP -g-MA-HNT/PP nano -composites . .......................... 8 Figure 1.3 . Schematic representation of serial sectioning (a) front view, and (b) top v iew . ......... 10 Figure 1.4 . Multiscale imaging of the HNT composite (a) slide generated through serial sectioning of HNT composite, (b) HNTs imaged with Auriga SEM and (c) individual HNT . ............................................................................................................................. 12 Figure 1.5 . (a) 2 -D SEM images of HNT polypropylene composite, (b) 3 -D reconstruction of the RVE based on serial sectioning . .............................................................................. 13 Figure 1.6 . Basic steps in the realization algorithm (OF=object ive function; MC =Monte Carlo) . . 14 Figure 1.7 . Schematic of nucleation and growth of the c ells during realization process . ............. 16 Figure 1.8 . Two -point correlation functions . .................................................................................. 16 Figure 1.9 . (a) The mesh after its final simplification for real RVE, (b) inclusions™ mesh in real inclusions™ mesh in statistical RVE . ............................................................................................................................... 17 Figure 1.10 . (a) SEM image of thermoplastic/xGnP composite, (b) TPCFs of thermoplastic/xGnP in three directions, (c) Realization of isotropically dispersed xGnP in polymer, (d) TPCFs of isotropically dispersed xGnP in polymer in three directions, (e) preferred oriented xGnP in polymer with constant radius=average radius of xGnP in SEM image, (f) TPCF curves of real mictrostructre and reconstructed one, (g) Extracting PDF of xGnP diameters from SEM images of nan ocomposite, (h) incorporating size distribution into construction cod e and making more accurate RVE . .................................................... 23 Figure 2.1 . Schematic of (a) microcomposite, (b) exfoliated, and (c) intercalated polymer Œclay nanocomposite morphologies (Bourbigot, VanderHart et al. 2003) ........................... 32 Figure 2.2 . (a) Single platelet, (b) Bulk powder, (c) Edge view and (d) Typical aggromeration (Drzal, Fu kushima et al. 2006) . ..................................................................................... 33 Figure 2.3 . Nano -Composite Processing (High Density Polyethylene/GNP) . ................................. 36 Figure 2.4 . Processing of high density p olyethylene/GNP composite . .......................................... 37 Figure 2.5 . SEM images of xGnP+HDPE . ........................................................................................ 39 Figure 2.6 . Gradient in xGnP distribution across the injected molded polymer nanocomposite .. 40 x Figure 2.7 . A typical load Œdisplacement curve (top) and (bottom) the deformation pattern of an elastic Œplastic material during and after indentation (bottom) (Oliver and Pharr 1992) . .......................................................................................................................... 42 Figure 2.8 . (a) xGnP/HDPE 3vf% flexural injected molded specimen, (b) nanocomposite specimen fixed in epoxy fixture, (c) na noindents and (d) SEM images ........................ 43 Figure 2.9 . Load -displacement curve in the middle of the specimen . ........................................... 44 Figure 2.10 . Load -displacement curv e near the edge . ..................................................................... 44 Figure 2.11 . Experimental setup for tensile test and DIC (digital image correlation) .. .................... 46 Figure 2.12 . Uniaxial tensile test on HDPE in two different strain rates . ......................................... 47 Figure 2.13 . Uniaxial tensile test on HDPE and HDPE/GNP, crosshead speed=0.8 mm/min . .......... 47 Figure 2.14 . (a) ASTM D638 dogbone specimen, (b), (c) proposed shear specimen (units in ch) .. 49 Figure 2.15 . Applied force versus elongation for uniaxial tensile samples . ..................................... 50 Figure 3.1 . 3D representation of the spherical coordinates of a randomly selected point . .......... 59 Figure 3.2 . Examples of 3D models of nanocomposites with different aspect ratios (AR), a) VF=1%, AR=100, b) VF=1%, AR=10 . .............................................................................. 59 Figure 3.3 . Typical traction -separation law for modeling cohesive failure . .................................. 61 Figure 3.4 . Comparison of traction Œseparation response in opening and sliding separation modes (Awasthi, Lagoudas et al. 200 9. ....................................................................... 65 Figure 3.5 . Averaged stress Œstrain curves for RVEs with random distribution and orientation of GNP; Particle volume fractions VF = 0.5, 1.5, 2.5%; Aspect ratio D/t = 40; perfectly bonded GNP . ................................................................................................................ 68 Figure 3.6 . Averaged stress Œstrain curves for RVEs with random distribution and orientation of GNP; Particle volume fractions VF = 1%; Aspect ratio D/t =10, 50, 100; perfectly bonded GNP . ................................................................................................................ 69 Figure 3.7 . Comparing GNP/HDPE (perfect bonding with different VF and constant AR) and GNP/HDPE (cohesive bonding w ith different VF and constant AR) ............................. 70 Figure 3.8 . Left: GNP/HDPE (perfectly bonded with different VF and constant AR); Right: GNP/HDPE (cohesively bonded with different VF and constant AR) ........................... 71 Figure 3.9 . Comparing GNP/HDPE (perfectly bonded with different AR and constant VF) and GNP/HDPE (cohesively bonded with different AR and constant VF) ........................... 72 xi Figure 3.10 . Left: GNP/HDPE (perfectly bonded with different AR and constant VF); Right: GNP/HDPE (cohesively bonded with different AR and constant VF) . .......................... 73 Figure 3.11 . Study showing the effect of AR (weak bonding) . ......................................................... 75 Figure 3.12 . Study showing the effect of VF (weak bonding) . ......................................................... 75 Figure 3.13 . Damage sequence of HDPE/GNP with AR=100, VF=1%, Left: weak bonding, Right: strong bonding . ............................................................................................................ 76 Figure 4.1 . Summary of pot entials used for CNT, Graphene and polyethylene . ........................... 87 Figure 4.2 . Generating PE chains using a self -avoiding walk algorithm . ........................................ 91 Figure 4.3 . (a) Nano RVE for CNT/PE, (b) Nano RVE for Gr/PE and (c) nano RVE for xGnP/PE . ..... 92 Figure 4.4 . Polymer chains around the fillers in equilibrium state . ............................................... 93 Figure 4.5 . Concentration profiles of polymer around Graphene sheet . ....................................... 94 Figure 4.6 . Displacement of individual carbon atoms of Graphene during the pull -out . ............. 95 Figure 4.7 . Increments of axial force applied to the n anotube and Graphene over time. Velocity of the center of mass of a carbon NT and Graphene during the simulated pullout .. .. 96 Figure 4.8 . Displacement of individual carbon atoms of Graphene and CNT during the pull -out process .......................................................................................................................... 96 Figure 4.9 . Axial tensile responses of Graphene+PE (1,5 and 10 Graphene layers) . ................... 97 Figure 5.1 . Void nucleation rate . .................................................................................................. 105 Figure 5.2 . Void evolution during uniaxial tension . ...................................................................... 106 Figure 5.3 . Void evolution during pure shear . .............................................................................. 107 Figure 5.4 . Void shearing (Weck, Wilkinson et al. 2006) . ............................................................. 108 Figure 5.5 . A schematic drawing illustrates the void shear mechanism. (Xue 2008) . .................. 108 Figure 5.6 . Stress -strain response (uniaxial tensile test) head speed=0.8 mm/min . ................... 114 Figure 5.7 . Plastic strain histories as a function of Poisson™s ratio compared with undamaged material (isotropic material) . ...................................................................................... 115 Figure 5.8 . GTN model, isotropic material (Von Mises) and experimental uniaxial tests . ........... 116 xii Figure 5.9 . Limit of the porosity ff at fracture as a function of stress triaxiality . ......................... 118 Figure 5.10 . Limit of the porosity ff at fracture for the balanced biaxial mode, effect of q1 . ....... 119 Figure 5.11 . Limit of the porosity ff at fracture for the balanced biaxial mode, effect of q1 . ....... 120 Figure 5.12 . uniaxial tensile tests, GTN model versus experimen tal results . ................................. 121 Figure 5.13 . Mesh of shear specimen . ........................................................................................... 123 Figure 5.14 . Load -displacement curve for shear specimen,experiment versus simulation . .......... 123 Figure 5.15 . Void evolution during loading in the modified specimen . ......................................... 124 Figure 5.16 . RVE model with 1% weight fraction of GNP, Finite element meshes of (a) RVE (b) GNP . ............................................................................................................................ 125 Figure 5.17 . Multipoint constraint boundary condition . ................................................................ 127 Figure 5.18 . Damage evolution in GNP/HDPE composite. (SDV2: void density) . ........................... 128 Figure 5.19 . Overall engineering stress -strain for GNP/HDPE composite versus HDPE . ............... 129 xiii KEY TO SYMBOLS Elastic modulus E Possion™s ratio Shear modulus G Volume fraction Stiffness tensor C Length of the RVE L Displacement along directions x,y,z u,v,w Displacement variables Stress variables Stain variables Volume Reaction force F Mass matrix M Force vector F Lame™s constants Density Triaxiality Effective Mises stress q Hydrostatic pressure p Material parameters in GTN model q1, q2 ,q3 Void volume fraction f xiv Stretch variables x y z GTN Gurson -Tvergaard -Needleman model GTN MD Molecular Dynamics MD Two -point correlation function TPCF Work done by the traction nG,sG,tG Effective traction at damage initiation eff T Maximum value of the effective separation max m Aspect ratio AR Yield function Accumulated effective strain The effective porosity Von Mises effective stress M ()YfMf1 Chapter 1 3-D microstructure reconstruction of polymer nano -composite using FIB ŒSEM and statistical correlation function Application to halloysite clay nanotubes polypropylene composite 1.1 Introduction and literature review In recent years, polymer nano -composites (PNCs) have increasingly gained more attention due to their improved mechanical, barrier, thermal, optical, electrical and biodegradable properties (Kojima, Usuki et al. 1993 , Messersmith and Giannelis 1994 , Sinha Ray and Okamoto 2003) in comparison with the conventional micro -composites or pristine polymer. With a modest addition of nanoparticles (usually less than 5 wt.% (Kojima, Usuki et al. 1993 )), PNCs offer a wide range of improvements in moduli (Okada and Usuki 1995 , Giannelis 1996 , Schmidt, Shah et al. 2002 ), strength, heat resistance (Yano, Usuki et al. 1997 ), biodegradability (Bharadwaj 2001, Sinha Ray and Okamoto 2003 ), as well as decrease in gas permeability (Messersmith and Giannelis 1994 , Gilman 1999 , Schmidt, Shah et al. 2002 , Schmidt and Malwitz 2003 ) and flammability (Haile 1992 , Messersmith and Giannelis 1994 , Bharadwaj 2001 , Schmidt, Shah et al. 2002 , Afaghi Khatibi and Mortazavi 2008 ). Although PNCs offer enormous opportunities to design novel material systems, development of an effective numerical modeling approach to predict their properties based on their complex multi -phase and multiscale structure is still at an early stage. Different experimental and simulation approaches are used to measure/calculate the thermomechanical properties in nanoscale. Molecular Dynamics (MD) is becoming a powerful computational tool for the simulation of matter at the molecular scale (Torquato 1997 , Nemat -Nasser and Hori 1999 , Komanduri, Chandrasekaran et al. 2001 , Pham and Torquato 2003 , 2 Mortazavi 2009 ). To estimate macro level or bulk properties of nano -composites, multisca le homogenization approaches are utilized based on continuum mechanics principles. The homogenization techniques can be categorized into the following six classes: statistical methods such as weak and strong -contrast (Dumont, Ladevèze et al. 1987 , Nemat -Nasser and Hori 1999 ), inclusion -based methods such as self -consistent or Mori ŒTanaka (Hori and Munasighe 1999 ), numerical methods such as finite element analysis and asymptotic methods (Affdl and Kardos 1976), variational/energy based methods such as Hashin ŒShtrikman bounds (Sheng, Boyce et al. 2004), and empirical/semi -empirical methods such as Halpin ŒTsai and classical upper and lower bounds (Voigt ŒReuss) (Zhu and Narh 2004 ). Finite element method for continuum mechanics has also been successfully applied to the integrated representative volume elements (RVEs) with a nanometric secondary phase (Hsueh, Fuller et al. 1999 , Fertig III and Garnich 2004 , Dong, Bhattacharyya et al. 2008 ). Regardle ss of the method to reproduce microstructures of the nano -composite, either with well aligned RVE (Hsueh, Fuller et al. 1999 , Dong, Bhattacharyya et al. 2008 ), or randomly distributed RVE with Monte Carlo scheme (Fertig III and Garnich 2004 ), the final RVE cannot entirely represent the actual complex and highly heterogeneous nano -composite structures. Dong et al. (Cannillo, Leonelli et al. 2002 ) developed a framework to incorporate the microstructural images such as SEM micrographs into 2 -D finite element modeling, the so called object -oriented finite element (OOF) technique. In their work, they combined data from the real microstructures such as particle size, shape, spatial position, and orientation distribution with fundamental material parameters inclu ding elastic modulus, Poisson™s ratio, and coefficient of thermal expansion (CTE) of the constitutive phases to understand the overall material behavior. The OOF however, is limited to elasticity and thermal conductivity calculations in two -dimensional 3 microstructures (Cannillo, Pellacani et al. 2003 , Wang, Kulkarni et al. 2003 , Velichko, Holzapfel et al. 2007 , Ray 2010 ). Several experimental and theoretical techniques such as X -ray computed tomography (CT), and focused ion beam/scanning electron microscopy (FIB/SEM) are being exploited to obtain three -dimensional microstructure as the input RVE for FEM software (Faessel, Delisée et al. 2005 , Baniassadi, Laachachi et al. 2011 , Deng, Liu et al. 2012 ). We have also used this method in the current work to obtain microstructure information of nano - composite. However, due to the high cost of sample preparation processes, simulation methods are often more desirable for the reconstruction of heterogeneous mic rostructures. Statistical reconstruction of heterogeneous materials using statistical correlation functions can be used as an alternative tool to reconstruct heterogeneous materials. Two -point correlation function (TPCF) is the simplest statistical corre lation functions that can convey some information about dispersion and distribution of inclusions in heteroge neous materials. Characterizati on and reconstruction of materials can be performed using different orders of statistical correlation functions (N -Point Correlation functions) (Dumont, Ladevèze et al. 1987 , Nemat -Nasser and Hori 1999 , Torquato 2002 , Fullwood, Adams et al. 2008 , Baniassadi, Garmestani et al. 2011 , Hamedani, Baniassadi et al. 2011 , Baniassadi, Ahzi et al. 2012 , Baniassadi, Mortazavi et al. 2012 ). In multiphase het erogeneous m aterials, the one point correlation function represents the volume fraction of different phases and does not describe any structural information (Baniassadi, Garmestani et al. 2011). TPCFs are the simplest well -known class of statistical descriptors that can be used to describe, in addition to volume fraction, the morphology and distribution of the heterogeneous materials (Baniassadi, Ahzi et al. 2012 ). Even higher order correlation functions must be used (Baniassadi, Ahzi et al. 2011 ) in order to 4 increase the precision of the statistical continuum approach. TPCFs can be measured from different techniques such a s microscopy (SEM or TEM) (Hamedani, Baniassadi et al. 2011 ), small X -ray scattering (Baniassadi, Laachachi et al. 2011 ), or Monte Carlo simulations (Baniassadi, Laachachi et al. 2011 ). Recently, Deng et al. (Den g, Liu et al. 2012 ) presente d a statistical work based on 2 -D realization of the microstructure obtained from SEM images of carbon black particle fillers disp ersed in synthetic natural rub ber. Their statistical approach is based on TPCF and two point cluster function using annealing technique. In this study, unlike previous studies, the 3 -D reconstruction of the microstructure of polypropylene nano -composites with 10 wt.% (7.2 vf.%) HNT fillers was ac hieved using (a) a three dimens ional (3 -D) morphology -based R VE, and (b) an RVE of nano -com posite constructed using statistical TPCFs. To our knowledge, this is the first time a statistical reconstructi on method based on two -dimensional microstructu re information of nano -composit e has been used to approximate the three -dimensional microstructure. Secondly, finite element analysis was used to deform the RVEs under tension and shear deformations to measure the effective stiffness tensor of the HNT polymer composite. Finally, the nu merically predicted results obtained from both RVEs were compared against the measured experimental data, in order to assess the feasibility of the statistical approach in predicting the various properties of anisotropic nano -composites. 1.2 Statistical correlation function The one -point correlation function of the phase q is defined by the probability of occurrence of random points in this phase. Mathematically, the one -point correlation function can be expressed as (Torquato 2002 ): 5 qiq11 Total NS(r) N (1) Where, N total is a large number of random points which are selected in the non -Eigen microstructure and N i represent the number of events which occur in phase q. TPCFs are determined based on the probability of occurrence of the head and tail of each vector in a particular phase. TPCFs in each direction are measured by selecting a number of random vectors such as r of the length rr within the microstructure and then calculating the probability of positioning the sta rting and ending points of the vector, i.e., T1 and T2, in desired phases (for better illustration see Figure 1.1). Depending on the state of the starting and ending points of r, for a m -phase medium there will be 2m probabilities, mathematically expressed as 2112 ,ij ij ijNNPrrrrTT N (2) Where ij N is the number of vectors with the beginning in phase (i, i), and the end in phase (j, j) with ,1,2,..., ijm . 1r and 2r stand for position vectors, respectively, 1T and 2T. It should be pointed out that by definition the TPCFs of anisotropic medium depend on the orientation of vector r. All TPCFs of a m -phase medium are not interdependent. For instance, in a three -phase composite there exist nine TP CFs consisting of 11P, 12P, 13P, 21P, 22P, 23P, 31P, 32P and 33P. Due to normality conditions the following equalities are immediately satisfied (Torquato 2002 ). 6 1,3 ijijPrv (3) 1,31,3 1ijijPr (4) i is volume fraction. With TPCFs being symmetric (in non -FGM microstructures), i.e. ijji PP, there remains only three independent TPCFs in the above example of a three -phase medium (Baniassadi, Ahzi et al. 2011 ) . In general, for a m -phase medium where 3m, there are 212mm independent TPCFs. For the HNT polymer composite, there exist exactly two states, phase -1 (po lymer matrix), and phase -2 (HNT particles). Therefore, only one of the four TPCFs is independent. Figure 1.1 . Four different configurations of vectors for two -point correlation function. In this work, the Monte Carlo simulation has been conducted to measure TPCFs of SEM image. In this approach, TPCFs are estimated by assigning a number of random vectors within 7 digitized SEM images, and examining the number fraction of the sets (vectors) w hich satisfy the desirable types of correlation states (Baniassadi, Laachachi et al. 2011 ). 1.3 Materials and synthesis The homopolymer polypropylene (PP) Pro -fax 6301 (LyondellBasell) and modifier (2% polypropylene -graft maleicanhydride, PP -g-MA) were first dried under vacuum overnight at 80 C. The molten PP pellets and the modifier were mixed using a co -rotating twin -scre w extruder (DSMXplore) , at 190 C for 3 min and extruded and chopped into pellets. The pellets were compounded with different weight fractions (5wt.%, 10 wt . % and 20 wt . %) of HNT (NaturalNanoInc.), at 190 C for 3min with a rota tion speed of 100 rpm. HNT polypropylene nanocomposite was processed without modifier as well. The HNT weight fraction used in this paper was 10 wt. % HNT/PP. 1.4 Mechanical properties The evaluation of mechanical properties of nano -composite and the host polymer were carried out using a UTS mechanical testing system and properties were measured based on ASTM standard dogbone specimens . To check the reproducibility of the experi - mental data and to ensure their consistency , 5 specimens were tested for each formulation . The tensile modu lus of the 10 wt.% HNT/PP and host polymer have been measured at 1.8 ± 0.03 GPa and 1.3 ± 0.04 GPa respectively. The effect of HNT loading on the tensile modulus of nano -composite with and without modifier is depicted in Fig ure 1.2 . This plot shows that by adding polypropylene graft maleicanhydride as a modifier , tensile modulus will be enhanced by 17%. Modifier enhances the bonding between HNT and PP. 8 Fig ure 1.2 . Tensile modulus of HNT/PP and PP -g-MA-HNT/PP nano -composites. 1.5 Serial sectioning of the nano -composite using FIB ŒSEM Simultaneous sectioning and imaging of the nano -composite (10 wt. % HNT+PP) was performed using a dual column focused ion beam (FIB) -SEM (CarlZeiss Auriga CrossBeam). Serial sectioning involved the removal of a known volume of the material by the ion beam followed by an incremental analysis with the electron beam. Because the sputtered material may redeposit onto the surface under analysis, significant in situ sample preparation was required. To begin, a trapezoid was milled int o the composite such that the shorter face was in a position to be imaged by the electron beam. The wider end of the trapezoid allowed for an unobstructed view of the analysis face. Two wings were on either side of the short face, such that after milling a shape similar to Fig ure 1.3 was observed. The wings were used as channels for sputtered 9 material to redeposit away from the surface of interest. A large beam (30kV, 20 nA) was used to excavate the bulk of the material and a smaller beam (30 kV, 4nA) was u sed to square the edges. The trenches were milled to a depth of 20 m. Water vapor was leaked into the chamber above the sample to assist the etching. A polished face was created by milling with a fine current beam (30kV, 1nA) to a depth of 20 m. A volume was then established in the software (SmartSEM, Carl Zeiss) with a width and height larger than the viewing area. A milling current of 1nA was used again. A schematic of the serial sectioning is shown in Fig ure 1.3 and the real images recorded during the FIB procedure are presented in Fig ure 1. 4. The width of each slice was 50 mm, therefore 50 nm of the nano - composite would be milled away with the ion beam followed by an image capture with the electron beam. The i mage contrast was turned slightly higher than what would normally be used to acquire a good image to accentuate the HNT from the matrix and aid in the reconstruction. Around sixty to one hundred slices were taken per sample, a process that took 2 Œ3 h. A se ries of 2D images representing slices or cross sections of the RVE is generated through FIB ŒSEM cutting. The advantage of using serial sectioning is to obtain a series of slides with the same reference point, allowing an automated 3D reconstruction techniq ue to be applied. 10 Fig ure 1. 3. Schematic representation of serial sectioning ( top ) front view, and (b ottom ) top view. 11 1.6 3-D Reconstruction using VCAT® software The sixty serial sectioning bitmap files obtained through serial sectioning are imported into VCAT software (released by V -CAD Program, RIKEN, Japan), the 3 -D nano -composite is represented with gray levels between the range of colors 0 Œ255 according to the image binarizat ion mode 8 -bit HSV (hue, saturation, value) color map (Fig ure 1.5 ). By choosing a color threshold of (0,0, 100), the im age part representation gives the best approximation of the dimensions of the cluster of HNTs inside the matrix. A Mask Property is associated with the matrix that will be a color value between 0 and 255. A noise reduction filter is applied in order to smooth the surface on which the Mask Property is applied. In this study, each of the two phases (matrix and filler) identified in the material was given a unique ID in order to distinguish between the phases inside the nano -composite. The resulting 3 -D representative volume element (RVE), shown in Fig ure 1.5 , possesses the most realistic features (size, shape, and distribution) of the actual nano -composite, most su itable for the calculation of its material properties . VCAT software offers a relatively simple function for mesh generation, refinement and simplification, volume data storage and data transfer. The Mask Property option in VCAT allows the user to assign m aterial properties to each phase of the material prior to the finite element analysis to predict the mechanical properties of the nano -composite. Inclusion volume fraction of the reconstructed RVE was calculated using Mask Property software to be 7.2 vf.%. This value is very close to the experimental data with the equivalent weight fraction. 12 (a) (b) (c) Figure 1. 4. Multiscale imaging of the HNT composite (a) slide generated through serial sectioning of HNT composite, (b) HNTs imaged with Auriga SEM and (c) individual HNT (www.naturalnano.com) 13 Figure 1. 5. (a) 2 -D SEM images of HNT polypropylene composite, (b) 3-D reconstruction of the RVE based on serial sectioning. 1.7 3-D Reconstruction of nano -composite using statistical correlation functions In this study, SEM images were used to me asure TPCFs for the polypropyle ne matrix filled with HNTs nano -composite. Also, a previously developed Monte Carlo methodology was used as a mean for the 3 -D reconstruction of the microstru cture of heteroge neous nano -composite, based on TPCFs (Baniassadi, Garmestani et al. 2011 ). The salient feature of the reconstruction methodology is the ability to realize the 3 -D microstructure from 2 -D SEM images. 14 1.7.1 Reconstruction procedure TPCFs can be exploited as a primary statistical description of heterogeneous materials to reconstruct three -dimensional microstructures. Each statistically reconstructed microstructure is called a fiRealizationfl. In the Monte Carlo reconstruction methodology adopted here (Torquato 2002), three -dimensional (3 -D) reconstruction of the microstructure is performed based on two -point statistical functions from tw o-dimensional (2 -D) SEM images. In this method, using the Monte Carlo algorithm (see figure 6), different realizations are generated to find a reconstructed microstructure that provides the best fit to the initial TPCFs. Figure 1. 6. Basic steps in the r ealization algorithm (OF=objective function; MC =Monte Carlo). 15 In the realization algorithm, different phases of the heterogeneous medium are represented by different cells, which are allowed to grow. Figure 7 shows the nucleation and growth of the cells during the realization process within 2D slices of the 3 -D reconstructed microstructure. The growth of cells, however, is controlled via several optimization parameters related to rotation, shrinkage, translation, distribution, and growth rates of cells (Garmestani, Baniassadi et al. 2009). Hybrid stochastic algorithm, based on the colony and kinetic algorithms, is a powerful capability that can be used to simulate wide variety of virtual microstruct ures. The simulation algorithm involves repeated realizations and minimization of an objective function (OF) at the end of the preceding realization. The OF is defined based on TPCFs derived from each of the realizations and the real microstructure. This a lgorithm is presented based on three steps of the realization process; 1) Generation, 2) Distribution, and 3) Growth of cells. In the generation procedure, initial geometries are assigned to each phase. During the distribution step, these basic geometries (cells) are distributed in the simulation step. After distribution of basic cell geometries, the growth of cells starts using cellular automata approach. Realization steps are repeated until an adequately realistic microstructure is developed as compared s tatistically to the true microstructure (Baniassadi, Garmestani et al. 2011 ). In figure (8), TPCFs of reconstructed and real microstructure are compared, and the results show good agreements between them. 16 Figure 1. 7. Schematic of nucleation and growth of the cells during realization process. Figure 1. 8. Two -point correlation functions. 17 1.8 Finite element mesh using vcat2tets software The 3 -D reconstructed microstructures were exported to the mesh generation software called ‚vcat2tets™ developed by VCAD® to create an appropriate finite element (FE) mesh. Inclusions have complex shapes, vary in size and orientations, and their surface cu rvature varies from one surface element to another. The distance between inclusions is not uniform; some inclusions are very closely spaced, whereas others are not. To take into account the influence of these geometric details on the local stress and strai n distributions, a very fine mesh has been used to define surfaces of the inclusion, and in the regions between closely spaced inclusions. The mesh was then simplified to a target number of tetrahedral elements were achieved. For example, the original FE m esh contained 25,000,000 tetrahedral elements, however, by using the ‚SimpTets™ software developed by VCAD® the final mesh size was reduced to 1000,000 tetrahedral elements (Figure 1.9). To insure that the structure is represented accurately by the simple mesh, finer mesh sizes were also tested and the mesh size independency was verified. Figure 1. 9. (a) The mesh after its final simplification for real RVE, (b) inclusions™ mesh in real RVE with average inclusion volume of 0.025 m3, and (c) inclusions™ mesh in statistical RVE. The 3 -D RVE obtained based on statistical method was also meshed using vcat2tets software. Figure 8 shows the FE mesh for both models after the final simplification. 18 1.9 Finite element analysis of RVEs The FE model illustrated in Figu re 8 is a rectangular cube representing a RVE with 5.06 x 5.22 x 4.15 3m volume. The created FE mesh was used for numerical simulation and calculation of the mechanical response of the material under variou s loading conditions. Tensile modulus of the host polymer measured in section 2.2 has been used (1.3±0.04 GPa) for matrix. Tensile modulus of HNT has been measured numerically and experimentally by many researchers. For example, Guimaraes, L., et al. (Guimarães, Enyashin et al. 2010 ) calculated the HNT modulus using molecular dynamic simulation and recently Dong Lu et al. (Lu, Chen et al. 2011) have measured the young™s modulus of halloysite nanotube using transmission electron with a bending stage. They have shown that the inner and outer diameters of the hall oysite nanotube are the key parameters to determine the tensile modulus of this nanoparticle. HNT which has been used in this study has the inner, outer diameters and length of 20nm and 85nm and 1000nm respectively. Based on the information available in t he work by Dong Lu et al. we have estimated the tensile modulus of HNT 140 GPa. Since the primary goal of this paper was to validate the viability of 3 -D reconstruction of the real microstructure using FIB/SEM and statistical methods, we used the simplest material model to characterize the nano -composite. It was therefore assumed that both polypropylene matrix and HNT would behave as linear elastic materials with perfect interfacial bonding between the two constituents. This idealization of the material be havior will be relaxed in the future, by assuming rate dependent materials as well as including a third interfacial phase between the matrix and the inclusion. 19 The FE simulations were carried out using ABAQUS TM commercial finite element software (ABAQUS In c., Providence, RI). A total of 6 simulations of the mechanical tests (3 tensile tests, and 3 shear tests) were carried out using this RVE. The objective was to find the modulus of elasticity and Poisson™s ratios of the RVEs in the X, Y and Z directions. The generalized form of Hooke's law, as shown in equation 5, was employed in order to perform the calculations. :S (5) In Eq. (5), the tensor S is the compliance tensor, and and are strain and stress tensors, respectively. One percent (1%) deformation strain was applied to each RVE. Using the stress and strain distributions for each test, the volumetric average value of stresses and strains were calculated using equations (6) and (7), in which N is the total number of elements. Volumetric stress: (6) Volumetric strain: (7) The compliance tensor ( S) was then constructed, with results shown in Eq. (8). Comparing the constructed compliance tensor for the RVE with the compliance tensor for an isotropic material, as shown in Eq. (9), this HNT composite was concluded to be an isotropic material. As a result, the diagonal terms, S 11 , S22 , and S 33 , give the in verse of the elastic modulus. The elastic modulus in X, Y, and Z directions were found to be, Exx=1.76 GPa, Eyy=1.53 GPa and Ezz=1.93 GPa, as shown in Table 1. NixxxiVolume iStress Stress 1)()(NixxxiVolume iStrain Strain 1)()(20 7.10059.0005.00017.0002.0007.0006.047.1004.0024.0016.0001.0006.00033.068.1014.0002.0004.00007.00163.0006.0516.015.0141.00025.000757.0004.0136.0653.0150.0002.00025.00003.0134.0162.0568.0S GPa -1 (8) 123123332211123123332211)1(2000000)1(2000000)1(20000001000100011E (9) The same procedure was also carried out with the statistical RVE and the results are summarized in Table 1.1. The average module of elasticity w as calculated to be 1.74 GPa and 1.64 GPa for real and statistical RVE, which are within 5% and 8.8 % of the exp erimental value of 1.8 +0.03 GPa. The difference between the estimated value using numerical approach and experimental value is related to the interphase region. Unlike microsized inclusions, interphase has a large influence on the overall properties of the nanocomposite. The perturbation of the polymer chains near nanoparticle creates a constrained region around the nanoparticle. (Kojima, Y.; Usuki, A.; Kawasumi, M.; Okada, A.; Fukushima,Y.; Kurauchi, T.; Kamigaito, O. J Mater Res 1993, 8, 1185.). This inte rphase region poses the same length scale as that of the nanoparticle but its properties is different from the host matrix. It has been shown by many researchers that taking into account the effect of interphase in evaluating the tensile modulus of nanocom posite will give us higher value for the modulus. 21 Table 1.1 . Properties of Halloysite clay nanotube and polypropylene (Lu, Chen et al. 2011 ) (top) and, elastic constants measured using real RVE and statistical RVE (bottom). Halloysite Density 2.5 g/cc Elastic Modulus 140 GPa Poisson Ratio 0.4 Polypropylene Density 0.9 g/cc Elastic Modulus 1.3+ 0.04 GPa Poisson Ratio 0.3 Real RVE Statistical RVE Exx 1.76 1.64 Elastic Modulus Eyy 1.54 1.64 (GPa) Ezz 1.93 1.64 xy 0.33 0.33 Poisson's Ratio yz 0.34 0.33 zx 0.33 0.34 1.10 3D Reconstruction of xGnP/Polymer using statistical correlation functions In this section, the s tatistical method developed for Halloysite nanotube has been used to construct the microstructure of anisotropic xGnP polymer composite, where only image of one slice of the actual microstructure has been used to obtain TPCF curves. Since the nanocomposite microstructure is not isotropic, TPCF curves in three different directions need to be matched with the corresponding curves of the real microstructure (see Figure 1.10.a, 1.10. b). The goal is to find a microstruc ture that has the same TPCF curves in these directions. An isotropically dispersed nanocomposite has been reconstructed using Monte Carlo method (Figure 1.10.c) and 22 TPCF curves in thre e directions are shown in F igure 1.10.d. TPCF curves in all three direct ions are the same, as expected due to the isotropy of the model. GNP nanoplat elets have been reoriented by controlling three Euler Angles to obtain the same TPCF curves of the real microstructure. The average GNP diameter obtained from SEM images has been used for GNP (see Figure 1.10.e) . Figure 1.10.f shows the best matching curves obtained using this approach. The probability distribution function of the GNP diameter has been obtained using multiple SEM images of nanocomposites and this function has been used in the reconstruction code (see Figure 1.10.g) . So instead of using constant diameter for GNP (average diameter), GNP diameter has a distribution obtained from e xperiment. Figure 1.10.h shows a good fit between TPCF curves of reconstructed and real microstructure. This microstru ctured inspired RVE can be used for mechanical, thermal or electrical characterizations of any platelet nanocomposite. In Chapter 3 inter facial debonding of GNP polymer nanocomposite will be studied using RVE obtained from this approach. 23 Figure 1. 10. (a) SEM image of thermoplastic/xGnP composite, (b) TPCFs of thermoplastic/xGnP in three directions, (c) Realization of isotropically dispersed xGnP in polymer, (d) TPCFs of isotropically dispersed xGnP in polymer in three directions, (e) preferred oriented xGnP in polymer with constant radius=average radius of xGnP in SEM image, (f) TPCF curves of real mictrostructre and reconstructed one, (g) Extracting PDF of xGnP diameters from SEM images of nanocomposite, (h) incorporating size distribution into construction code and making more accurate RVE 00.020.04 0.06 0.080.10.12 0.1402004006008001000TPCF r 0 /4 /8 Realization in MATLAB 00.020.040.06 0.080.10.12 0.14050100TPCF r 0 Degree 22.5 Degree 45 Degree Constant Radius= Average Radius 00.02 0.04 0.06 0.08 0.1 0.12 0.14 010002000TPCF r (a) (b) (c) (d) (e) (f) 24 Figure 1.10 (cont™d) 1.11 Conclusion and future work The primary goal of this research was to assess the feasibility of two separate methods for the 3-D reconstruction of the actual microstructure of nano -composites. To investigate the accuracy of these two methods, elastic mechanical properties of the HNT p olypropylene composite were obtained based on two models; RVE of the real material™s micro/nanostructures constructed using morphological images that were captured with scanning electron microscopy (SEM), and RVE obtained based on statistical methods. The modeling techniques were compared with experimental data. This showed that the RVE of the nano -composite obtained based on TPCFs can partially approximate real microstructure because TPCF cannot capture all statistical 00.00020.00040.00060.00080.001001000200030004000Diameter (nm)PDF00.02 0.04 0.06 0.08 0.1 0.12 0.14 0200 400 600 800 1000 TPCF r (g) (h) 25 information of the nano -composite (Gommes, Jiao et al. 2012 ). But in future we will add other additional statistical descriptors to have better approximation for reconstructi on. Fortunately, the developed reconstruction procedure is very flexible to add other statistical descriptors such as lineal path (Manwart, Torquato et al. 2000 ). The 3 -D reconstruction using statistical approach has several advantages over the 3 -D reconstruction based on FIB ŒSEM. FIB ŒSEM is time intensive, expensive, and not applicable to every material. For instance, our attempts to perform FIB ŒSEM on carbon bas ed nano -composites, e.g., exfoliated graphite nanoplatelets (xGnP) polymer composite, were unsuccessful due to the poor contrast in SEM images caused by the small difference in molecular weight of xGnP and the host polymer. We are working to rectify this p roblem by expanding the statistical method developed for Halloysite nanotube for the anisotropic xGnP polymer composite, where only one slice or image of the actual microstructure will be sufficient to reconstruct the nano -composite. Although in this paper , we investigated a nano -composite with only one volume fraction, the TPCF of a nano -compo site with a higher or lower volume fractions can be estimated using the method developed by PIs in (Ghazavizadeh, Soltani et al. 2012 ). This will greatly facilitate future optimization of nano -composites with optimum volume frac tion, and reduces the need for obtaining SEM images. In this paper, it was also assumed that polypropylene matrix and clusters of HNT have perfect interfacial bonding. However, in the case of nano -composites, interphase properties can have a significant im pact on material properties. Therefore, in the future, following the improved statistical method presented in this paper, an RVE comprised of three phases; xGnP, interphase, and the host polymer, will be reconstructed in order to account for the effect of interphase properties on mechanical properties of nano -composites. 26 REFERENCES 27 REFERENCES Afaghi Khatibi, A. and B. Mortazavi (2008). "A Study on the Nanoindentation Behaviour of Single Crystal Silicon Using Hybrid MD -FE Method." Advanced Materials Research 32: 3. Affdl, J. and J. Kardos (1976). "The Halpin -Tsai equations: a review." Polymer Engineering & Science 16(5): 344 -352. Baniassadi, M., et al. (2011). "New approximate solution for N -point correlation functions for heterogeneous materials." Journal of the Mechanics and Physics of Solids . Baniassadi, M., et al. (2012). "New approxima te solution for N -point correlation functions for heterogeneous materials." Journal of the Mechanics and Physics of Solids 60(1): 104 -119. Baniassadi, M., et al. (2011). "Three -phase solid oxide fuel cell anode microstructure realization using two -point correlation functions." Acta materialia 59(1): 30 -43. Baniassadi, M., et al. (2011). "Mechanical and thermal behavior of nanoclay based polymer nanocomposites using statistical homogenization approach." Composites Science and Technology 71(16): 1930 -1935. Baniassadi, M., et al. (2011). "Statistical continuum theory for the effective conductivity of carbon nanotubes filled polymer composites." Thermochimica Acta 520(1-2): 33 -37. Baniassadi, M., et al. (2012). "Three -dimensional reconstruction and homogenization of heterogeneous materials using statistical correlation functions and FEM." Computational Materials Science 51(1): 372 -379. Bharadwaj, R. K. (2001). "Modeling the barrier properties of polymer -layered silicate nanocomposites." Macromolecules 34(26): 9189 -9192. Cannillo, V., et al. (2002). "Numerical models for thermal residual stresses in Al2O3 platelets/borosilicate glass ma trix composites." Materials Science and Engineering A 323(1-2): 246-250. Cannillo, V., et al. (2003). "Numerical modelling of the fracture behaviour of a glass matrix composite reinforced with alumina platelets." Composites Part A: Applied Science and Manufacturing 34(1): 43 -51. Deng, H., et al. (2012). "Utilizing Real and Statistically Reconstructed Microstructures for the Viscoelastic Modeling of Polymer Nanocomposites." Composites Science and Technology . Dong, Y., et al. (2008). "Experimental chara cterisation and object -oriented finite element modelling of polypropylene/organoclay nanocomposites." Composites Science and Technology 68(14): 2864 -2875. 28 Dumont, J., et al. (1987). "Damage mechanics for 3 -D composites." Composite structures 8(2): 119-141. Faessel, M., et al. (2005). "3D Modelling of random cellulosic fibrous networks based on X -ray tomography and image analysis." Composites Science and Technology 65(13): 1931 -1940. Fertig III, R. S. and M. R. Garnich (2004). "Influence of constituent properties and microstructural parameters on the tensile modulus of a polymer/clay nanocomposite." Composites Science and Technology 64(16): 2577 -2588. Fullwood, D. T., et al. (2008). "A strong contrast homogenization formulation for multi -phase anisotro pic materials." Journal of the Mechanics and Physics of Solids 56(6): 2287 -2297. Garmestani, H., et al. (2009). "Semi -inverse Monte Carlo reconstruction of two -phase heterogeneous material using two -point functions." International Journal of Theoretical and Applied Multiscale Mechanics 1(2): 134 -149. Ghazavizadeh, A., et al. (2012). "Composition of two -point correlation functions of subcomposites in heterogeneous materials." Mechanics of materials . Giannelis, E. P. (1996). "Polymer layered silicate nanocomposites." Advanced materials 8(1): 29-35. Gilman, J. W. (1999). "Flammability and thermal stability studies of polymer layered -silicate (clay) nanocomposites1." Applied Clay Science 15(1-2): 31 -49. Gommes, C. J., et al. (2012). "Microstructural degeneracy associated with a two -point correlation function and its information content." Physical Review E 85(5): 051140. Guimarães, L., et al. (2010). "Structural, electronic, and mechanical propertie s of single -walled halloysite nanotube models." The Journal of Physical Chemistry C 114(26): 11358 -11363. Haile, J. M. (1992). Molecular dynamics simulation : elementary methods . New York, Wiley. Hamedani, H. A., et al. (2011). "Microstructure, propert y and processing relation in gradient porous cathode of solid oxide fuel cells using statistical continuum mechanics." Journal of Power Sources 196(15): 6325 -6331. Hori, M. and S. Munasighe (1999). "Generalized Hashin ŒShtrikman variational principle for boundary -value problem of linear and non -linear heterogeneous body." Mechanics of materials 31(7): 471 -486. Hsueh, C., et al. (1999). "Analytical and numerical analyses for two -dimensional stress transfer." Materials Science and Engineering: A 268 (1): 1 -7. 29 Kojima, Y., et al. (1993). "Mechanical properties of nylon 6 -clay hybrid." Journal of Materials Research(USA) 8(5): 1185 -1189. Komanduri, R., et al. (2001). "Molecular dynamics (MD) simulation of uniaxial tension of some single -crystal cubic metals at nanolevel." International Journal of Mechanical Sciences 43(10): 2237-2260. Lu, D., et al. (2011). "Direct Measurements of the Young's Modulus of a Single Halloysite Nanotube Using a Transmission Electron Microscope with a Bending Stage." Journal of Nanoscience and Nanotechnology 11(9): 7789 -7793. Manwart, C., et al. (2000). "Stochastic reconstruction of sandstones." Physical Review E 62(1): 893. Messersmith, P. B. and E. P. Giannelis (1994). "Synthesis and characterization of layered silicate -epoxy nanocomposites." Chemistry of materials 6(10): 1719 -1725. Mortazavi, B., Khatibi, A. A., Politis, C. (2009). "Molecular Dynamics Investigation of Loading Rate Effects on Mechanical -Failure Behaviour of FCC Metals." journal of computational and theoretical nanoscience 6: 644 -652. Nemat -Nasser, S. and M. Hori (19 99). Micromechanics: overall properties of heterogeneous materials , Elsevier Amsterdam. Okada, A. and A. Usuki (1995). "The chemistry of polymer -clay hybrids." Materials Science and Engineering: C 3(2): 109 -115. Pham, D. and S. Torquato (2003). "Strong -contrast expansions and approximations for the effective conductivity of isotropic multiphase composites." Journal of applied physics 94(10): 6591-6602. Ray, S. S. (2010). "A new possibility for microstructural investigation of clay -based polymer nanoco mposite by focused ion beam tomography." Polymer 51(17): 3966 -3970. Schmidt, D., et al. (2002). "New advances in polymer/layered silicate nanocomposites." Current Opinion in Solid State and Materials Science 6(3): 205 -212. Schmidt, G. and M. M. Malwitz (2003). "Properties of polymer -nanoparticle composites." Current opinion in colloid & interface science 8(1): 103 -108. Sheng, N., et al. (2004). "Multiscale micromechanical modeling of polymer/clay nanocomposites and the effective clay particle." Polyme r 45(2): 487 -506. Sinha Ray, S. and M. Okamoto (2003). "Polymer/layered silicate nanocomposites: a review from preparation to processing." Progress in polymer science 28(11): 1539 -1641. 30 Torquato, S. (1997). "Effective stiffness tensor of composite medi a--I. Exact series expansions." Journal of the Mechanics and Physics of Solids 45(9): 1421 -1448. Torquato, S. (2002). Random heterogeneous materials : microstructure and macroscopic properties . New York, Springer. Velichko, A., et al. (2007). "3D chara cterization of graphite morphologies in cast iron." Advanced Engineering Materials 9(1-2): 39 -45. Wang, Z., et al. (2003). "Effects of pores and interfaces on effective properties of plasma sprayed zirconia coatings." Acta materialia 51(18): 5319 -5334. Yano, K., et al. (1997). "Synthesis and properties of polyimide -clay hybrid films." Journal of Polymer Science Part A: Polymer Chemistry 35(11): 2289 -2294. Zhu, L. and K. Narh (2004). "Numerical simulation of the tensile modulus of nanoclay -filled polyme r composites." Journal of Polymer Science Part B: Polymer Physics 42(12): 2391 -2406. 31 Chapter 2 Processing and testing of nanocomposites 2.1 Processing of nanocomposites The final properties of nanocomposites mainly depend on the processing techniques (Manias, Touny et al. 2001 , Hussain, Hojjati et al. 2006 , Jang and Zhamu 2008 ). In order to take advantage of the outstanding properties of nano -fillers , the nano fillers must be well -dispersed in the polymer matrix. The higher the degree of exfoliation, the better the nanofiller can transfer the load, leading to a higher stiffness. However, nanofillers tend to agglomerate or aggregate because of the small size and large surface areas so it is very difficult to uniformly distribute nanofillers. Thus , a processing technique that results in a better dispersion of nano -particles would lead to better properties in the resulting nanocomposite. There are three different morphologies of nanocomposites based on the degree of exfoliation of the nanofiller, microcomposites, intercalated na nocomposites and exfoliated nanocomposites (Figure 2.1). In microcomposites, there is no penetration of the polymer into the nano platelet. In an intercalated nanocomposites, polymer penetrates into the nanoplatelet structure . In exfoliated nanocomposites, the individual nanoplatelets are dispersed as single platelets into the polymer matrix. 32 Figure 2.1. Schematic of (a) microcomposite, (b) exfoliated, and (c) intercalated polymer Œclay nanocomposite morphologies (Bourbigot, VanderHart et al. 2003 ) 2.2 Processing of polymer/ graphene nanocomposites Graphene nanoplatelet (GNP) is a new class of carbon nanoparticles with multifunctional properties. GNP has fiplateletfl morphology, meaning they are very thin but with wide aspect ratio. This wide dimension and platelet morpho logy makes them an excellent barrier to liquid and gases, while their pure graphitic composition makes them outstanding electrical and thermal conductors (Fukushima, Drzal et al. 2006 ). Advantages of GNP in mechanical reinforcement over existing carbon fillers have also been addressed in numerous literature (Cipriano, Kota et al. 2008, Kim and Macosko 2009 , Steurer, Wissert et al. (2009) ). Figure 2.2 shows xGnP or exfoliated graphite nano -platelet used in this study. Layered Clay Monomer Polymer Polymer Microcomposite Exfoliated nanocomposite Intercalated nanocomposite 33 Figure 2.2. (a) Single platelet, (b) Bulk powder, (c) Edge view and (d) Typical aggromeration (Drzal, Fukushima et al. 2006 ) Recently, polymer/graphene nanocomposites have drawn great attention. Gra phene can be dispersed in several polymer matrices; these nanocomposites exhibit significantly improved mechanical and thermal properties due to the dispersion of low weight fraction loadings of nanometer -sized layered graphene with high aspect ratio and h igh strength in the polymer matrix. GNP has a similar layered structure as nano -clay, the processing techniques for polymer/nano -clay can be used for GNP nanocomposites such as melt intercalation , solution mixing , and in -situ polymerization. In the melt intercalation method, a solid mixture of nanoplatelets and polymer is melted up to melting point of the polymer so polymer chains can penetrate, exfoliate and intercalate the nanoplatelets (Ray and Okamoto 2003 , Deshmane, Yuan 2 µm 1 1 1 µm 1 µm 5 nm (a) (b) (c) (d) 34 et al. 2007 ). In the solution mixing method, polymer first is dissolv ed in a solvent and then nanoplatelets will be mixed in the solution (Ray and Okamoto 2003 ). In In-situ polymerization method , first nanoplatetes are immersed in a monomer solution and then m onomer penetrate s into the nanoplatetes . Monomer will be polymerized inside the interlayer nanoplatelets galleries and will be exfoliated. This method can be used for limited polymers such a s nylon 6, poly (methyl methacrylate, PMMA), epoxy, and phenolic resin (Yeh, Liou et al. 2004 , Pappas, Patel et al. 2005 ). In-situ polymerization and solution mixing results in a b etter mechanical, thermal and electrical properties in nanocomposites due to a better dispersion and stronger interactions between nano -reinforcements and the polymer matrix. However, the conventional compounding method of melt -mixing followed by injection molding or compression molding is still considered to be the major processing method used for manufacturing the thermoplastics in the industry because of its design flexibility, low cost and labor (Ibeh 2014 ). 2.3 Processing of high density polyethylene -exfoliated graphene nanoplatelet nanocomposites In this study, melt extrusion follow ed by injection molding was employed to manufacture high density polyethylene (HDPE) Œ exfoliated graphene nanoplatelet (GnP) nanocomposties. HDPE pellets with the trade name Marlex® HXM 50100 (Density 0.95 g/cm3, Flow index 10.0 g/10 min,MW ~230,000) were obtained from Chevron Phillips Che mical Company. Paraffin wax (max C30, Density 0.92 g/cm3, MW~ 500) with the melting point of 55 º C was purchased from Sigma ŒAldrich. GNP nanoplatelets were obtained from XG Science, Inc. The diameter of individual GNP platelet is around 15 µm (GNP -15) and the thickness is less than 10 µm. The physical properties of GNP -15 is detailed in Table 2.1. 35 Melt extrusion of HDPE/GNP nanocomposites was carried out in a DSM Micro 15 cc Compounder, (Vertical, co-rotating, twin -screws microextruder) operating at 220 ºC for 5 min at a screw speed of 100 rpm. The composite melt was then transferred to a Daca Micro injector with the Tbarrel = 220 ºC and Tmold = 90 ºC. The injection pressure applied for the injection molding of flexural coupons was at 0.6 MPa. The melt extrusion and injection molding systems are shown in Fig ure 2. 3. Table 2.1. Geometrical and surface characteristics of GNP -15 (Jiang and Drzal 2010 ) Material Length (µm) Diameter (µm) Aspect ratio Surface area (m2/g) Intrinsic density (g/cm3) GnP -15 <0.01 (platelet thickness) 15 (platelet diameter) ~1500 ~40 2.1 Large scale extrusion was also performed using a Leistritz extruder with 25 mm co-rotating twin screws and a module multilayer slit die provided by Premier Dies Corporation, which yields a 16 layer structure . The extruder and die were set to 200ºC with a s crew rotation speed of 20 rpm and yielding a melt pressure from 2000 -2300 psi depending on the xGnP concentration. Upon exiting the die, the extrudate was cooled with a 3 roll chill stack set to 100 ºC, and finally collected on a spool at a rate of 2 feet per minute Figure 2. 4. 36 xGnP/HDPE Extruder Injection Molding Paraffin Wax Xylenes Mixing Wax HDPE xGnP Sonication for required energy(100 Watt) and time (20 min) Figure 2.3. Nano -Composite Processing (High Density Polyethylene /GNP) 37 Figure 2.4. Processing of h igh density p olyethylene /GNP composite HDPE GNP Cold water 38 2.4 Scanning electron microscopy Scanning electron microscopy (SEM) was used to determine particle size and distribution and to examine fracture surfaces. The SEM consists of an electron gun producing a source of electrons at an energy range of 1 -40keV. Electron lenses reduce the diameter of the electron beam and place a small focused beam on the specimen. The electron beam interacts with the near - an image. The smaller the beam size, the better the resolution of the image. Operating the SE M requires fine tuning to optimize picture quality with resolution. SEM is run under a vacuum to minimize beam interactions with gas molecules which would retard resolution. Non -conductive specimens, such as most polymers, often suffer from variations in s urface potential which introduce astigmatism, instabilities, and false x -ray signals. Charging, a condition during which charge accumulates on the surface of a non -conducting specimen causing excessive brightness, often occurs making it difficult to obtain quality images. Sputter coating non -conductive samples with a fine gold layer is often required to avoid these issues. The primary goal of using SEM was to determine particle size and particle dispersion. SEM images were taken of the dogbone fracture surf aces. The fracture surfaces were sputter coated with gold atoms prior to imaging to avoid charging, but that did not negatively affect the image quality. Figure 2.5, exhibit s example images showing the dispersion of particles. Particle size was measured fr om such images and has been used for calculating the probability density function for GNPs in Chapter 1. 39 Figure 2 .5. SEM image s of xGnP+HDPE 10 µm 10 µm 10 µm 10 µm 20 µm 10 µm 20 µm 100 µm 40 2.5 Nanoindentation The spatial variation of the superstructure developed along and across the flow direction of injection -molded HDPE/GNP composite is presented using SEM imaging. T here is a significant gradient in mechanical, thermal and electrical properties across the nanocomposite specimen due to the different level of nanoparticle dispersion in different region. Shear stress near mold surfaces is higher than the one in the middle of the s pecimen so this affects the level of dispersion and exfoliation of the nanoplatelets . This has also been observed by other researchers for particulated polymer composite (Yalcin and Cakmak 2004 , Konishi and Cakmak 2005 ). SEM images have been taken from several regions of specimen and results have been shown in Figure 2.6. As it is shown in this figure , xGnP platelets are well dispersed near mold walls and are less dispersed in the middle since shear stress near mold wall is higher. Figure 2. 6. Gradient in xGnP distribution across the injected molded polymer nanocomposite Nanoindentation was used to determine the elastic modulus of nanocomposite in the middle and near the edges of the sample . The material undergoes both e lastic and plastic deformation under the indenture in nanoindentation , but upon unloading elastic deformation will be 41 recovered . Figure 2.6 shows a typical load -displacement curve and the deformation pattern of an elastic Œplastic material during and after indentation . Oliver and Pharr 1992 related t he initial unloading contact stiffness to the reduced modulus, Er, according to the following equa tion . 2rSAE A is the projected contact area, and S is the experimentally measured stiffness, S, which is determined as a slope ( S=dP/dh, P is the load, h is the displacement) of the upper portion of the unloading curve (Figure 2.7). The reduced modulus is given by 22(1) 1(1) iriEEE In which E and are the Young™s modulus and Poisson™s ratio for the specimen, and E i and i are the same parameters for the Berkovich indenter (E i=1141 GPa and i =0.07). The hardness is calculated using th is equation (Oliver and Pharr 1992 ) max PHA Pmax is the peak indentation load. The unloading curve follows power low ()mfPBhh Where B and m are empirically evaluated fitting parameters. According to Oliver and Pharr , the area function for a perfect Berkovich indenter is given by 2()24.5 cc Ahh where h c is the vertical distance along which contact is made 42 max max cPhh S and is a constant that depends on the indenter geometry; for a Berkovich indenter =0.75. Figure 2. 7. A typical load Œdisplacement curve (top) and (b ottom ) the deformation pattern of an elastic Œplastic material during and after indentation (botto m) (Oliver and Pharr 1992 ). Load, P hc for =1 hc for =0.72 Loading Unloading Possible range for hc hmax Surface profile after removing load Indenter Initial surface Surface profile under load 43 A Nano Indenter XP (MTS Systems Co.) was used to carry out the indentation measurements with the DCM techniques. The theoretical displacement resolution and load resolution are 0.0002nm and 1 nN, respectively, and the maximum loading is 10 mN. The Nano Indenter is operated by the MTS TestWorks software enabling instrument control and data acquisition and processing. A Berkovich indenter with a radius of 150 nm was used for all the experiments. The contact area function was fit for contact depths ranging from 5 to 500nm using curve fit weighting at shallow depths. The DCM was used at a dyna mic frequency of 75 Hz . DCM measurements were carried out at an amplitude of displacement oscillations of 2 nm, a strain rate of 0.05 s -1 , and a drift rate of <0.05 nms -1. Each specimen was tested in 25 indentations arranged in a matrix of 5×5 indentation s with separation distances of 25µ m (see figure 2. 8). Figure 2. 8. (a) xGnP/HDPE 3vf % flexural injected molded specimen, (b) nanocomposite specimen fixed in epoxy fixture, (c) nanoindents and (d) SEM images a. Nano -composite 3% vf b. Nano -composite fixed in epoxy c. Nanoindents d. SEM 12 mm 3 mm 50 m 50 m 44 The results of the nanoindenta tion have been shown in figures 2. 9 and 2. 10. Modulus of elasticity has been calculated based on these curv es an d the results are in table 2.2 and 2.3.As it was expected elastic modulus is lower in the middle compare to elastic modulus near edges due to the higher volume fraction of GNP near the edges. Figure 2.9. Load -displacement curve in the middle of the specimen Table 2.2. Elastic modulus measured in the middle of nanocomposite Mean (Elastic modulus GPa) 1.973 Std. Dev. 0.602 0510 1520 25 303502000 4000 6000 Load (mN) Displacement into surface (nm) 45 Figure 2.10.Lload -displacement curve near the edge Table 2.3. Elastic modulus measured near the edge of nanocomposite Mean (Elastic modulus GPa) 1.77 Std. Dev. 0.398 2.6 Mechanical tests A mechanical testing machine ( United Calibration Corp SFM 20 ) was used to measure the tensile properties, according to ASTM D638 standard. System control and data analysis were performed using the Datum software. A 1000 lbs load cell was used for measuring the tensile strength and tensile modulus of the injection mo lded tensile coupons. Digital image correlation (DIC) technique was utilized to precisely calculate the strain history. Images recorded from a digital camera positioned to monitor deformation in the gauge area of each specimen are used to compute the strain. Figure 2.11 shows the experimental set -up. Crosshead speeds of 0.2, 0.8 and 2 mm. min-1 were used for testing the plastics. All results presented are the average values of five measurements. Tensile test has been performed in two different test speeds to show the polymer 051015 20 25 30350200040006000Load (mN) Displacement into surface (nm) 46 strain rate dependency. As it is shown in Figure 2.12 thermoplastic shows a stiffer response under higher strain rate. Tensile tests also have been performed on HDPE/GNP samp les and the result is shown in Figure 2.13. There is a significant improvement in stress -strain response of the final product by adding only 1wt% GNP in to polymer. Both stiffness and ultimate tensile strength have been improved. Figure 2.11. Experimental setup for tensile test and DIC (digital image correlation) . 47 Figure 2.12. Uniaxial tensile test on HDPE in two different strain rates Figure 2.13. Uniaxial tensile test on HDPE and HDPE/GNP , crosshead speed=0.8 mm/min 01020 304000.10.20.30.4True stress (MPa) Strain Neat Polymer,strain rate effect 2mm/min 0.2mm/min 0510 152025 3035404500.1 0.2 0.3 0.4 0.5 True stress MPa True strain HDPE HDPE+1wt% GNP HDPE+2.2 wt% GNP 48 2.7 Shear test Shear tests are of great importance since they are widely used to investigate the shear modulus and shear strength of materials. In this study, a specially designed specimen is used to perform the shear stress tests in a tensile testing machine to get the results in terms of forces and extension, crack initiation etc. A modified shear test specimen combining the simple uniaxial test with a zone of interest is tested that almost provides the pure shear . The dogbone tensile specimens processed using the big i njected molded machine are machined and t ensile and proposed modified shear test specimens are shown in Figure 2.14. Shear specimen is loaded in uniaxial tension to find the correlation data points between force and extension (Figure 2.1 5). Test speed is 0.2 mm/min with data acquisition rate of 10 Hz. This load displacement data will be used in calibration of damage model for HDPE in Chapter 5. The whitening due to void formation in HDPE has been shown in Figure 2.14 c. 49 Fig ure 2.1 4. (a) ASTM D638 dogbone specimen, (b), (c) proposed shear specimen (units inch) Whitening due to void (a) (b) (c) 50 Figure 2.15. Applied force versus elongation for uniaxial tensile samples 0306090120 150 180 0246Load (N) Extension (mm) Force -disp. for shear test 51 REFERENCES 52 REFERENCES Bourbigot, S., et al. (2003). "Investigation of nanodispersion in polystyrene Œmontmorillonite nanocomposites by solid -state NMR." Journal of Polymer Science Part B: Polymer Physics 41(24): 3188-3213. Cipriano, B. H., et al. (2008). "Conductivity enhancement of carbon nanotube and nanofiber -based polymer nanocomposites by melt annealing." Polymer 49(22): 4846 -4851. Deshmane, C., et al. (2007). "High strength Œtoughness combination of melt intercalated nanoclay -reinforced thermoplastic olefins." Materials Science and Engineering: A 460: 277 -287. Drzal, L., et al. (2006). "Exfoliated graphite nano platelets (xGnP): a carbon nanotube alternative for modifying the properties of polymers and c omposites." XG Sciences . Fukushima, H., et al. (2006). "Thermal conductivity of exfoliated graphite nanocomposites." Journal of thermal analysis and calorimetry 85(1): 235 -238. Hussain, F., et al. (2006). "Review article: polymer -matrix nanocomposites, processing, manufacturing, and application: an overview." Journal of composite materials 40(17): 1511 -1575. Ibeh, C. C. (2014). Thermoplastic Materials: Properties, Manufacturing Methods, and Applications , CRC Press. Jang, B. Z. and A. Zhamu (2008). " Processing of nanographene platelets (NGPs) and NGP nanocomposites: a review." Journal of Materials Science 43(15): 5092 -5101. Jiang, X. and L. T. Drzal (2010). "Multifunctional high density polyethylene nanocomposites produced by incorporation of exfoli ated graphite nanoplatelets 1: Morphology and mechanical properties." Polymer Composites 31(6): 1091 -1098. Kim, H. and C. W. Macosko (2009). "Processing -property relationships of polycarbonate/graphene composites." Polymer 50(15): 3797 -3809. Konishi, Y. and M. Cakmak (2005). "Structural hierarchy developed in injection molding of nylon 6/clay/carbon black nanocomposites." Polymer 46(13): 4811 -4826. Manias, E., et al. (2001). "Polypropylene/montmorillonite nanocomposites. Review of the synthe tic routes and materials properties." Chemistry of Materials 13(10): 3516 -3523. Oliver, W. C. and G. M. Pharr (1992). "An improved technique for determining hardness and elastic modulus using load and displacement sensing indentation experiments." Journa l of materials research 7(06): 1564 -1583. 53 Pappas, J., et al. (2005). "Structure and properties of phenolic resin/nanoclay composites synthesized by in situ polymerization." Journal of Applied Polymer Science 95(5): 1169 -1174. Ray, S. S. and M. Okamoto (2003). "Polymer/layered silicate nanocomposites: a review from preparation to processing." Progress in polymer science 28(11): 1539 -1641. Steurer, P., et al. (2009). "Functionalized graphenes and thermoplastic nanocomposites based upon expanded graphite oxide." Macromolecular rapid communications 30(4-5): 316 -327. Yalcin, B. and M. Cakmak (2004). "Superstructural hierarchy developed in coupled high shear/high thermal gradient conditions of injection molding in nylon 6 nanocomposites." Polymer 45(8): 2691-2710. Yeh, J. M., et al. (2004). "Comparative studies of the properties of poly (methyl methacrylate) Œclay nanocomposite materials prepared by in situ emulsion polymerization and solution dispersion." Journal of Applied Polymer Science 94(5): 1936 -1946. 54 Chapter 3 An interfacial debonding -induced damage model for graphite nanoplatelet polymer composites 3.1 Introduction and literature survey Graphene nanoplatelet (GNP) is a new class of carbon nanoparticles with multifunctional properties. GNP has fiplateletfl morphology, meaning they are very thin but with wide aspect ratio. This wide dimension and platelet morphology makes them an excellent ba rrier to liquid and gases, while their pure graphitic composition makes them outstanding electrical and thermal conductors (Fukushima, Drzal et al. 2006 ). Advantages of GNP in mechanical reinforcement over existing carbon fillers have also been addressed in numerous literature (Cipriano, Kota et al. 2008, Kim and Macosko 2009 , Ste urer, Wissert et al. (2009) ). Recently, polymer/graphene nanocomposites have drawn great attention. Graphene can be dispersed in several polymer matrices; these nanocomposites exhibit significantly impr oved mechanical and thermal properties due to the dispersion of low weight fraction loadings of nanometer -sized layered graphene with high aspect ratios and high strengths in the polymer matrix. In situ tensile tests show damage initiates in polymer nanoco mposites mainly by interfacial debonding. In an experimental study Rafiee et al. 2010 showed that functionalized graphene sheets are significantly effective at improving the fracture energy, fracture toughness, strength, stiffness, and fatigue resistance of epoxy polymers at remarkably lower nanofiller volume 55 fractions in comparison to carbon nanotubes (CNTs) and nanoclay additives. This can be related to their high specific surface area, two -dimensional geometry, and strong nanofiller -matrix adhesion. In another experimental work, Achaby et al. (El Achaby and Qaiss 2012 ) compared the effect of carbon based nanofiller on the tensile properties of the resulting nanocomposite. Their results showed that at the same filler volume fraction, the GNP performs better than the CNT, due to the higher specific surface area, larger aspect ratio and nanoscale 2 -D flat surface of the GNP. These properties of the GNP result in an enhanced mechanical interlocking with the polymer chains, and an enlarged interphase zone at the nanofiller polymer interface. Conducting experiment at the nanoscale, in order to understand the micromechanics of nanocomposites is difficult, if not sometime impossible. Therefore, computational and analytical methods must be used to study the me chanics of nanocomposites. A deep understanding of the damage and fracture mechanisms of nanocomposites is crucial for structural design and practical applications. Although damage mechanisms of traditional composites have been widely studied in the litera ture (McCartney 1987 , Orifici, Herszberg et al. 2008, Zhang, Yang et al. 2010 , Kushch, Shmegera et al. 2011 , Kushch, Shmegera et al. 2011 , Bheemreddy, Chandrashekhara et al. 2013 ), there are few studies (Zerda and Lesser 2001 , Wang, Chen et al. 2005 , Yoshimura and Okabe 2011 , Bheemreddy, Chandrashekhara et al. 2013 ) that report on the damage and fracture mechanisms of nanoco mposites. Theoretical and numerical predictions of the effective mechanical properties of fiber or particle reinforced nanocomposites are usually made under the assumption of high interfacial strength or perfect bonding. However, the interface behavior can significa ntly affect the 56 mechanical properties of nanocomposites, therefore an assumption of strong or perfect bonding would be inappropriate for these types of composites. Some researchers investigated nanofillers debonding in nanocomposites using analytical appr oaches. For example, Williams (Williams 2010 ) and Lauke (Lauke 2008 ) analyzed the energy dissipation phenomena b y considering, besides particle debonding, voiding and subsequent yielding of the polymer. Salviot et al. (Salviato, Zappalorto et al. 2013 ) developed a hierarchical analytical multi -scale model to assess the fracture toughness improvements due to the debonding of nanoparticles and the plastic yielding of nanovoids. Recently, they have also proposed a multiscale analytical model to quantify th e toughness improvement due to the shear banding around nanoparticles (Zappalorto, Salviato et al. 2012).Their results showed that nanocomposite toughening is strongly affected by the size of nanoparticles and by surface treatments. In an analytical study, Zhang et al. (Zhang, Zhao et al. 2013 ) investigated the da mage for a 3D model of a single tube of nano carbon and polymer by assuming a cohesive model for the interface. Their results showed that the peak value of the macroscopic stress -strain curve is defined by the strength of the cohesive interface. This means that the higher the cohesive strength is, the larger the value of the macroscopic peak strength will be. Pisano et al. (Pisano, Priolo et al. 2012 ) studied the effect of the gallery failure mechanism on the macroscopic behavior of intercalated epoxy Œclay nanocomposites using the 2D representative volume element concept and finite element simulations. This effect was studied for different clay content s, aspect ratios and orientations, and with different fracture properties assigned to the galleries. The main result was that the gallery failure is the main mechanism for 57 strength reduction in the intercalated nanocomposite. They also found that nanocompo sites™ strength decreases with increasing clay content, which is in agreement with available experimental results (Dorigato, Morandi et al. 2012 ). Needleman et al. (Needleman, Borders et al. 2010 ), in a finite element study, investigated the effects of interphase thickness and interfacial strength between carbon nanotube and polymer on the stiffness and strengt h of the nanocomposite. They modeled interface between the polymer and the CNT by a phenomenological cohesive relation. In this chapter, a hierarchical multiscale model is developed to study the damage initiation in GNP/high density polyethylene composites . The cohesive zone model (CZM) has been adopted to capture the nanofillers debonding. It is worth noting that choosing appropriate cohesive parameters is the most important part in the modeling of debonding in nanocomposites . Therefore, the information ab out interfacial properties of GNP and polymer has been obtained from molecular dynamics (MD) simulations. A representative volume element (RVE) composed of GNPs and polymer matrix was created to study the overall stress -strain response of the nanocomposite . The main goal of this research work was to perform a systematic computational study on the effects of nanofillers/polymer bonding conditions on the macroscopic response of GNP/polymer composites for different GNP volume fractions, aspect ratio, and inte rfacial strength. 58 3.2 Nanocomposite model 3.2.1 Representative volume element (RVE) A 3D representative volume element (RVE) was created for the nanocomposite consisting of GNP and polymer. The RVE was generated using an in -house developed C ++ algorithm. Implementation steps used for developing the RVE with the Monte Carlo methodology are defined below. Numerical simulations were carried out inside a cubic unit cell of constant side length of 1000 units (units may be equally interpreted as nm). GNPs were modeled as simple discs dispersed inside the RVE. The geometry of each GNP was modeled as two parallel circular plates separated by the thickness of the GNP. Each circular plate in the volume of the RVE was identified by a normal vector, a cent er, and a radius. To achieve a uniform random scatter of GNPs using the Monte Carlo method, the center of each GNP was selected randomly inside the sample RVE. Then, the associated normal vector was specified by means of random homogeneous functions, to pr oduce uniformly distributed random points on the surface of a sphere, following 2cos(21) vArcu (1) 3.1, and u, v are random variables belonging to [0, 1]. The normal vectors thus selected, guarantee a uniform random distribution of GNP orientations. For generating each GNP, the procedure of random selection of its center and normal direction was followed successively and then the next GNP was identic ally created. 59 Figure 3.1. 3D representation of the spherical coordinates of a randomly selected point The optimum size of the RVE for each volume fraction and aspect ratio was determined by increasing the volume of the RVE until the homogenized stress -strain values no longer changed significantly. Figure 3.2 shows RVEs of nanocomposites with different aspect ratios. Figure 3.2. Examples of 3D models of nanocomposites with different aspect ratios (AR), a) VF=1%, AR=100, b) VF=1%, AR=10 3.2.2 Cohesive zone model The behavior of GNPs and the matrix interface is represented by cohesive zone model (CZM) defined in terms of bilinear traction/ separation law (Alger 1997 ). This model is implemented in commercial finite element software ABAQUS 6.10. C ohesive behavior can be surface based or element based. Damage is defined as a material property for the cohesive element but as an interaction property for the cohesive surface, even though the constitutive (a) (b) 60 relations are almost the same for both models. T o define the cohesive law, three parameters were used here: stiffness, strength and the critical value of the energy release rate. In this study, the cohesive surface model was used to model the debonding in GNP/Polyethylene composite. This damage process requires accurate characterization of the interfacial or bond material. However, in the GNP/polymer nanocomposites, the interaction between GNP and polymer is difficult to determine through experimental measurements. Thus, we used the results of the molecu lar dynamics simulation (MD) by Awasthi et al. (Awasthi, Lagoudas et al. 2009 ) on the interfacial interaction between graphene and polyethylene. Damage modeling simulates the degradation and eventual failure of the bond between two cohesive surfaces. The failure mechanism co nsists of two ingredients: a damage initiation criterion and a damage evolution law. The initial response is assumed to be linear. Figure 3.3 shows a typical traction -separation response with a failure mechanism. However, it is important to recognize that damage in the surface -based cohesive behavior is an interaction property and not a material property. Concepts of strain and displacement are reinterpreted as contact separations; contact separations are the relative displacements between the nodes on the slave surface and their corresponding projection points on the master surface along the contact normal and shea r directions. Stresses are defined for the surface -based cohesive behavior as the cohesive forces acting along the contact normal and shear directions divided by the current area at each contact point. 61 Figure 3.3- Typical traction -separation law for modeling cohesive failure The constitutive relation between these traction stresses and separation is given by (Version 2011) nnnsnt nsssst ntsttt nn ss tttkkk TtkkkK tkkk (2) Here, t RnR is the traction stress in the normal direction, tRsR, tRtR are traction stresses in the first shear and second shear directions, respectively, K is the nominal stiffness matrix, is the separation in the normal direction, and , are separations in the first shear and second shear directions, respectively. If the normal and shear components are decoupled, Eq . (2) reduces to, T max sep No separation Complete separation Damage starts here. 62 000000nnssttnnsstttkTtkK tk (3) After completing the linear elastic traction -separation, damage will be started. There are different criteria for damage initiation. Considering cohesive parameters obtained from the molecular dynamics simulation, the maximum stress criterion was used in t his study (see Eq (4) ). Based on this criterion, damage is assumed to initiate when the maximum contact stress ratio reaches the value of one. This criterion can be represented as: maxmaxmax ,,1 nstnst tttMAX ttt (4) In the above equation, the ramp function 0.5nnn ttt defines the normal contact stressin pure normal mode, st and tt are shear contact stresses along the first and the second shear direction, respectively. Once the corresponding initiation criterion is reached damage starts to evolve. The damage evolution law describes the rate at which the cohesive stiffness is degraded. A scalar damage variable, D, represents the overall damage at the contact point. It initially has a value of 0, if damage evolution is modeled; D monotonically evolves from 0 to 1 upon further loading after the initiation of damage. The contact stress compo nents are affected by the damage according to 63 (1),0 nnnnDtt ttothewise (5) (1) sstDt (6) (1) tttDt (7) Where nt, st and tt are the contact stress components predicted by the elastic traction -separation behavior for the current separations without damage. The dependence of the fracture energy on the mixed -mode can be defined based on a power law fracture criterion. The power law criterion stat es that failure under mixed -mode conditions is governed by a power law interaction of the energies required to cause failure in the individual (normal and two shear) modes. It is given by 1nst CCCnst GGG GGG (8) With the mixed -mode, the fracture energy is equal to Cnst GGGG when the above condition is satisfied. In the above expression, the quantities nG, sG and tG refer to the work done by the traction and its conjugate separation in the normal, the first, and the second shear directions, respectively. The quantities of CnG, CsG, and CtG, which re fer to the critical fracture energies required to cause failure in the normal, the first, and the second shear directions should (Reeder 2006 , Dai and Mishnaevsky Jr 2013 ). For the 64 linear softening (see Figure 3.4), an evolution of the damage variable (Bheemreddy, Chandrashekhara et al. 201 3), D, reduces to : max0 max 0()()fmmm fmmm D (9) as the effective traction at damage initiation (defined eff Twith 2Cfmeff GTWhere to the maximum value of the effective separation attained during the refers max m below). .loading history 22effnst TTTT (10) 222mnst (11) Separation studies are conducted for both normal (traction) as well as sliding (shear) modes in the paper by Awasthi et al. (Awasthi, Lagoudas et al. 2009 ) and the cohesive zone parameters such as peak traction and energy of separation are evaluated for each mode. Traction -separation curves for normal debonding and sliding modes are shown in Figure 3.4. This figure shows that the maximum traction for normal mode is higher than the one for sliding mode while the separation point is higher in the shear mode. Cohesive parameters are listed in Table 1. 65 Figure 3.4. Comparison of traction Œseparation response in opening and sliding separation modes )Awasthi, Lagoudas et al. 2009 Table 3.1. Cohesive zone model parameters for opening and sliding modes Fracture mode Fracture energy P (mj/m P2P) Peak traction (MPa) Shear mode 331.650 108.276 Normal mode 246.525 170.616 3.3 Material behavior and FE simulation In this study, GNPs were assumed to behave as an isotropic linear -elastic solid with a modulus of E=1TPa , (Reddy, Rajendran et al. 2006 , Li, Dichiara et al. 2013 ). Normal Shear 66 An elastic -plastic model was used to model the high density polyethylene (HDPE) at low strain rate (Kwon and Jar 2008 ). The rate sensitivity was not considered for HDPE since the effect of the strain rate on the mechanical response of the polymer is only important at a high loading rate at which adiabatic deformation occurs (Dasari and Misra 2003 ). The required elastic -plastic model™s parameters were Young™s modulus, Poisson™s ratio, and a hardening curve (true tensile stress vs true tensile plastic strain). With the assumption of isotropic work hardening, the von Mises yield criterion was used. The von Mises criterion then relates the equivalent stress, e, to the yield stress in tension T by: eT (12) 222 122331 1()()() 2e (13) 1 2 3 are the principal stress components, and the tensile yield stress T is a material parameter with a minimum value, which denotes the limit of elastic behavior and the start of plastic deformation and will increase with tensile plastic strain. The assumed elastic properties for the linear isotropic polymer matrix were, E=11 00 MPa and .44. The equivalent von Mises stress for HDPE has been selected from [30] in a tabular form. A Python script was written to transfer the RVE based on the procedure explained in section 2.1 for ABAQUS/Explicit 6.10 commercial FE software . Three dimensional Six -node wedge element C3D6 and three -dimensional 4 -node linear tetrahedron elements C3D4 were used for GNPs and RVE, respectively. All RVEs were subjected to the uniaxial displacement loading on 67 one face while the opposite face was fixed. The computational time for the explicit FE simulation was reduced by mass scaling. In all simulations, the kinetic energy was less than 5% of the total strain energy indicating a quasi -static loading condition. Homogenized variables at the macroscopic scal e were obtained by volume averaging of variables in the RVE. The macroscopic stress and strain can be calculated from: 1mmmd (14) 1mmmd (15) Where and RmR RmR are the local (microscopic) stress and strain vectors and is the volume of the RVE. 3.4 Results and discussion 3.4.1 The effect of the GNP™s volume fraction and aspect r atio in perfectly bonded nanocomposite RVEs of nanocomposites with three different volume fractions (0.5, 1.5 and 2.5) have been created. To study the effect of volume fraction, the aspect ratio and the diameter of the inclusions were kept constant for all cases as; AR=40 and Dia.=10µm. As is shown in Figure 3.5, stiffness increases with an increase in the volume fraction. This is in agreement with the rule -of-mixture (Alger 1997 , Askeland, Fulay et al. 2011 ). 68 For the analysis of the effect of aspect ratio, RVEs of nanocomposites with three different aspect ratios (10, 50 and 100) have been created. The volume fraction and the diameter of the inclusions were kept constant for all cases as; VF=1% and Dia.=4 µm. The results show that the nanocomposites™ stiffn ess increases as the GNP™s aspect ratio increases (Figure 3.6). This is in agreement with the results obtained by Mortazavi et al. (Mortazavi, Baniassadi et al. 2013 ). These authors have compared the elastic modulus and thermal conductivity of two -phase random composites with different inclusion types and aspect ratios using finite elements and Mori ŒTanaka methods. They concluded that the elastic modulus and the thermal conductivity of nanocomposites increase by increasing the aspect ratio. Figure 3.5. Averaged stress Œstrain curves for RVEs with random distribution and orientation of GNP; Particle volume fractions VF = 0.5, 1.5, 2.5%; Aspect ratio D/t = 40; perfectly bonded GNP 0510 1520 2530 35404500.05 0.1 0.15 Stress(MPa) Strain AR=40,perfect bonding VF=.5 VF=1.5 VF=2.5 Pure matrix 04 8121600.004 0.008 0.012 Stress ( MPa) Strain 69 Figure 3.6. Averaged stress Œstrain curves for RVEs with random distributi on and orientation of GNP; Particle volume fractions VF = 1%; Aspect ratio D/t =10, 50, 100; perfectly bonded GNP 3.4.2 Comparing the effect of the GNP™s volume fraction and aspect ratio in perfectly bonded and cohesively bonded nanocomposites Here, we compare our predicted results for perfectly bonded and cohesively bonded composites. The cohesive zone model parameters are given in Table 1 and are based on the results of Awasthi et al. (see Figure 3.4). Our predicted results for the effects of volume fraction, with a constant aspect ratio of AR=40, are reported in Figure 3.7 and Figure 3.8. Figures 3.7a, 3.7b, and 3.7c show the effect of the volume fraction on the stress -strain responses of the perfectly bonded and cohesively bonded GNP/polymer composites in comparison to a pure polymer matrix. In all cases, when debonding starts to occur, the stress -strain curve for the damaged nanocomposite drops down to lower stress levels, and deviate from the perfectly bonded nanocomposite. This is in agree ment with the results of the recent 0102000.01 Stress ( MPa) Strain Perfect bonding 70 micromechanical model presented by Dastgerdi et al. (Nafar Dastgerdi, Marquis et al. 2014 ). Figure 3.8 also shows that as volume fraction increases, the difference between cohesively bonded and perfectly bonded responses will increase. This means that in nanocomposite with high volume fraction of GNPs, more inclusions will debond compared to th ose with low volume fraction. Figure 3.7. Comparing GNP/HDPE (perfect bonding with different VF and constant AR) and GNP/HDPE (cohesive bonding with different VF and constant AR) 0204000.1 Stress(MPa) Strain VF=0.5%,AR=40 pure matrix cohesive interface perfect bonding 0204000.1 Stress( MPa) Strain VF=1.5%,AR=40 Pure matrix cohesive interface perfect bonding 0204000.1 Stress( MPa) Strain VF=2.5%, AR=40 pure matrix cohesive interface perfect bonding 71 Figure 3.8. Left: GNP/HDPE (perfectly bonded with different VF and constant AR); Right: GNP/HDPE (cohesively bonded with different VF and constant AR) For the analysis of the effect of aspect ratio, we selected the composite of VF=1% and reported our predicted result s in Figure 3.9 and Figure 3.10 for different aspect ratios (AR=10, 50, 100). Figures 3.9a, 3.9b, and 3.9c show the effect of the aspect ratio on the stress -strain response of the perfectly bonded and cohesively bonded GNP/HDPE nanocomposites as they 0102030405000.05 0.1 0.15 Stress(MPa) Strain AR=40, perfect bonding VF=.5 VF=1.5 VF=2.5 Pure matrix 0102030 4000.05 0.1 0.15 Stress(MPa) Strain AR=40,Non -bonded VF=.5 VF=1.5 VF=2.5 Pure matrix 72 compa re with the response of a pure polymer matrix. In all cases, when the debonding starts to occur, the stress -strain curve for the damaged nanocomposite drop down to lower stress levels, and deviate from the perfectly bonded nanocomposites. Figure 3.9 also s hows that as AR increases, the difference between bonded and non -bonded response will increase. This is due to the high surface area in the higher aspect ratio which leads to the larger fracture surface during the debonding. Figure 3.10 shows that by incre asing the aspect ratio in the cohesive model, stress -strain response does not increase significantly as compared to the perfect bonding case. Figure 3.9. Comparing GNP/HDPE (perfectly bonded with different AR and constant VF) and GNP/HDPE (cohesively bonded with different AR and constant VF) 0204000.05 0.1 0.15 Stress (MPa) Strain VF=1%,AR=10 Perfect bonding cohesive interface pure matrix 0204000.1 0.2 Stress (MPa) Strain VF=1%,AR=50 cohesive interface perfect bonding pure matrix 0204000.1 0.2 Stress (MPa) Strain VF=1%, AR=100 Cohesive_interface perfect bonding pure matrix 73 Figure 3.10. Left: GNP/HDPE (perfectly bonded with different AR and constant VF); Right: GNP/HDPE (cohesively bonded with different AR and constant VF) 3.4.3 The effect of the GNP™s aspect ratio and volume fraction in weakly bonded nanocomposite To analyze the effect of w eak bonding, the cohesive zone parameters for weak bonding have been selected as shown in Table 2. The effect of the aspect ratio on the stress -strain response of the weakly bonded interface is plotted in Figure 3.11 for a fixed fillers volume fraction of 1% and compared with the response of the pure polymer matrix. This Figure indicates that: a) in the first stage of deformation, the stiffness of the nanocomposite increases as the aspect ratio of the nanofillers increase; and b) in the second part of defo rmation (after yield), the composite will have a lower flow stress compared to the host polymer. This indicates that debonding starts in the nonlinear region (high strains). Figure 3.11 also shows that the increase in the aspect ratio of the platelets resu lts in a lower flow stress of the composite. In fact, with a fixed volume fraction, nanocomposites with higher aspect ratio inclusions will have more defects compared to those 0204000.05 0.1 0.15 0.2 Stress(MPa) Strain Perfect bonding AR=100 AR=50 AR=10 pure matrix 01020304000.1 0.2 Stress(MPa ) Strain VF=1,Non -bonded AR=100 AR=50 AR=10 Pure matrix 74 with lower aspect ratio inclusions. Although by adding nanoplatelets to polymers , it is expected to have improved stiffness and toughness, the end result highly depends on the type of bonding generated between the GNP and the host polymer. To improve the interfacial bonding between nanofillers and the host polymer, different types of chemical modifiers have been investigated (Miyagawa and Drzal 2004 , Azeez, Rhee et al. 2013 ). The effect of the volume fraction of the fillers on the overall stress -strain response of the nanocomposite with weak interfacial bonding was shown in Figure 3.12 for a constant aspect ratio of AR=40. These results show that by increasing the volume fracti on of the nanofiller, the nanocomposite will have lower stress -strain response. This is due to the poor bonding between the filler and the matrix that causes more inclusion debonding for higher fillers volume fractions. Table 3.2. Cohesive parameters for weak bonding Fracture mode Fracture energy (mj/m 2) Peak traction (MPa) Shear mode 331.650 30 Normal mode 246.525 40 75 Figure 3.11. Study showing the effect of AR (weak bonding) Figure 3.12. Study showing the effect of VF (weak bonding) In order to clearly observe the damage sequence of HDPE/GNP nanocomposites from numerical predictions viewpoint, the state of damage at the interface during loading at strains of 0510 15 20 25 303500.05 0.1 0.15 Stress(MPa) Strain VF=1%, Weak bonding AR=10 AR=50 AR=100 pure matrix 01000.01 Stress(MPa) Strain VF=1 0510152025 303500.05 0.1 0.15 Stress(MPa) Strain AR=40, Weak bonding AR=10 VF=0.5% VF=1.5% VF=2.5% Pure matrix 76 6, 12 and 24 % i s shown in Figure 3.13. As one can observe, the interface damage variable D (called CSDMG) has reached a maximum value of 1.0 for many nodes present in the weakly bonded nanocomposites compared to the strongly bonded nanocomposites. Figure 3.13. Damage sequence of HDPE/GNP with AR=100, VF=1%, Left: weak bonding, Right: strong bonding e11 e11 e11 e11 e11 e11 77 3.5 Summary and conclusions In this chapter, a hierarchical multiscale model was developed to study the damage initiation in polymer/GNP nanocomposites. The cohesive zone model was used to model the polymer/nanofiller debonding. Results from atomic simulations of GNP pullout and debo nding tests were used for the cohesive zone model (CZM). Effects of volume fraction, aspect ratio, and fiber -matrix interfacial strength on the overall stress -strain response of the nanocomposite were investigated. Nanocomposites with perfectly bonded fill ers were also modeled for comparison with the nanocomposites with cohesive bonding. Results from studying the effect of fillers volume fraction (VF) and aspect ratio (AR) in perfectly bonded composites showed that both stiffness and toughness will increase with increasing VF and AR. The effect of the GNPs™ volume fraction on perfectly bonded and cohesively bonded composites was compared. As expected, the results showed that when the debonding starts to occur, the stress -strain curve for damaged nanocomposit e decreases and deviate from those for the perfectly bonded nanocomposite. Results also showed that as volume fraction increases, the difference between cohesively bonded and perfectly bonded responses will increase. This implies that in nanocomposites wi th higher volume fraction, more inclusions will debond ed compared to lower volume fraction. The effect of aspect ratio on the stress -strain response of the perfectly bonded and cohesively bonded GNP/HDPE nanocomposites has been studied. It was also shown that as AR increases, the difference between bonded and cohesively bonded response will increase. The cohesive model has also been tested for weakly bonded composites. The results showed by having weak bonding between inclusion and polymer the resulting c omposite will have a lower stress -strain response after the yield compared to the host polymer and that the increase in the aspect ratio of the platelets will lower the flow stress. In fact with fixed volume fraction, nanocomposites with 78 higher aspect rati o will introduce larger defect in the composite compared to nanocomposites with lower aspect ratio. Although by adding nanoplatelets into polymers, their stiffness and toughness will be improved provided good bonding exists between the nanofillers and the host polymer. Weak bonding between nanoplatelets and the polymer matrix may result in a composite with low stiffness and toughness. Hence, using different types of chemical modifiers for nanofillers to improve the interfacial bonding is necessary. We shoul d note that our simulation results show that the stiffness still increases with filler volume fraction in the case of weak bonding which indicates that debonding becomes noticeably important after the yield point. In this study it has been assumed that th e damage starts at the interface between nanofillers and the matrix. In chapter 5, matrix cracking following inclusion debonding will be modeled to capture a realistic damage behavior in nanocomposites. 79 REFERENCES 80 REFERENCES Alger, M. S. (1997). Polymer science dictionary, Springer. Askeland, D. R., et al. (2011). The Science and Engineering of Materials: Si Edition, CengageBrain. com. Awasthi, A. P., et al. (2009). "Modeling of graphene Œpolym er interfacial mechanical behavior using molecular dynamics." Modelling and Simulation in Materials Science and Engineering 17(1): 015002. Azeez, A. A., et al. (2013). "Epoxy clay nanocomposites Œprocessing, properties and applications: A review." Composites Part B: Engineering 45(1): 308 -320. Bheemreddy, V., et al. (2013). "Modeling of fiber pull -out in continuous fiber reinforced ceramic composites using finite element method and artificial neural networks." Computational Materials Science 79: 663 -673. Cipriano, B. H., et al. (2008). "Conductivity enhancement of carbon nanotube and nanofiber -based polymer nanocomposites by melt annealing. " Polymer 49(22): 4846 -4851. Dai, G. and L. Mishnaevsky Jr (2013). "Damage evolution in nanoclay -reinforced polymers: A three -dimensional computational study." Composites Science and Technology 74: 67 -77. Dasari, A. and R. Misra (2003). "On the strain rate sensitivity of high density polyethylene and polypropylenes." Materials Science and Engineering: A 358(1): 356 -371. Dorigato, A., et al. (2012). "Effect of nanoclay addition on the fiber/matrix adhesion in epoxy/glass composites." Journal of Composi te Materials 46(12): 1439 -1451. El Achaby, M. and A. Qaiss (2012). "Processing and properties of polyethylene reinforced by graphene nanosheets and carbon nanotubes." Materials & Design. Fukushima, H., et al. (2006). "Thermal conductivity of exfoliated graphite nanocomposites." Journal of thermal analysis and calorimetry 85(1): 235 -238. Kim, H. and C. W. Macosko (2009). "Processing -property relationships of polycarbonate/graphene composi tes." Polymer 50(15): 3797 -3809. Kushch, V., et al. (2011). "Numerical simulation of progressive debonding in fiber reinforced composite under transverse loading." International Journal of Engineering Science 49(1): 17 -29. 81 Kushch, V., et al. (2011). "E xplicit modeling the progressive interface damage in fibrous composite: Analytical vs. numerical approach." Composites Science and Technology 71(7): 989 -997. Kwon, H. and P. -Y. Jar (2008). "On the application of FEM to deformation of high -density polyeth ylene." International Journal of Solids and Structures 45(11): 3521 -3543. Lauke, B. (2008). "On the effect of particle size on fracture toughness of polymer composites." Composites Science and Technology 68(15): 3365 -3372. Li, W., et al. (2013). "Carbo n nanotube Œgraphene nanoplatelet hybrids as high -performance multifunctional reinforcements in epoxy composites." Composites Science and Technology 74: 221-227. McCartney, L. (1987). "Mechanics of matrix cracking in brittle -matrix fibre -reinforced compos ites." Proceedings of the Royal Society of London. A. Mathematical and Physical Sciences 409(1837): 329 -350. Miyagawa, H. and L. T. Drzal (2004). "The effect of chemical modification on the fracture toughness of montmorillonite clay/epoxy nanocomposites. " Journal of adhesion science and technology 18(13): 1571 -1588. Mortazavi, B., et al. (2013). "Modeling of two -phase random composite materials by finite element, Mori ŒTanaka and strong contrast methods." Composites Part B: Engineering 45(1): 1117-1125. Nafar Dastgerdi, J., et al. (2014). "Micromechanical modeling of nanocomposites considering debonding of reinforcements." Composites Science and Technology. Needleman, A., et al. (2010). "Effect of an interphase region on debonding of a CNT reinforced polymer composite." Composites Science and Technology 70(15): 2207 -2215. Orifici, A., et al. (2008). "Review of methodologies for composite material modelling incorporating failure." Composite structures 86(1): 194 -210. Pisano, C., et al. (2012). "Prediction of strength in intercalated epoxy Œclay nanocomposites< i> via finite element modelling." Computational Materials Science 55: 10 -16. Rafiee, M. A., et al. (2010). "Fracture and fatigue in graphene nanocomposites. " small 6(2): 179 - 183. Reddy, C., et al. (2006). "Equilibrium configuration and continuum elastic properties of finite sized graphene." Nanotechnology 17(3): 864. Reeder, J. R. (2006). "3D mixed -mode delamination fracture criteria Œan experimentalist™s perspective." 82 Salviato, M., et al. (2013). "Plastic shear bands and fracture toughness improvements of nanoparticle filled polymers: a multiscale analytical model." Composites Part A: Applied Science and Manufacturing. Steurer, P., et al. (2009). "Func tionalized graphenes and thermoplastic nanocomposites based upon expanded graphite oxide." Macromolecular rapid communications 30(4-5): 316 -327. Version, A. (2011). "6.11 Documentation." Dassault Syst emes Simulia Corp., Providence, RI, USA. Wang, K., et al. (2005). "Epoxy nanocomposites with highly exfoliated clay: mechanical properties and fracture mechanisms." Macromolecules 38(3): 788 -800. Williams, J. (2010). "Particle toughening of polymers by plastic void growth." Composites Science and Technol ogy 70(6): 885 -891. Yoshimura, A. and T. Okabe (2011). "Damage Growth Analysis in Particle -Reinforced Composite Using Cohesive Element." Advanced Composite Materials 20(6): 569 -583. Zappalorto, M., et al. (2012). "A multiscale model to describe nanocom posite fracture toughness enhancement by the plastic yielding of nanovoids." Composites Science and Technology. Zerda, A. S. and A. J. Lesser (2001). "Intercalated clay nanocomposites: morphology, mechanics, and fracture behavior." Journal of Polymer Sci ence Part B: Polymer Physics 39(11): 1137-1146. Zhang, B., et al. (2010). "A virtual experimental approach to estimate composite mechanical properties: Modeling with an explicit finite element method." Computational Materials Science 49(3): 645 -651. Zhang, Y., et al. (2013). "An analytical solution on interface debonding for large diameter carbon nanotube -reinforced composite with functionally graded variation interphase." Composite structures 104: 261 -269. 83 Chapter 4 Molecular dynamics simulation of polymer nanocompo sites 4.1 Introduction Due to the exceptional (mechanical, electrical and thermal) properties and its cost effectiveness, Graphene, a 2D counterpart of carbon nanotube (CNT), is an excellent candidate for nanocomposites. Graphene polymer nanocomposites have potential application in many technological fields such as electronics, packaging, sensors and membranes for gas separation. Recently, it has been revealed that the high mechanical performance of the polymer nanocomposites not only depends on the properties of the nanof iller, but also on the efficiency of load transfer from the polymer matrix to the nanofillers, which is determined by the interfacial bonding between them. Standard fiber pull -out tests could provide proof for strong interface between the graphene and the host matrix, however, these experiments are very expensive and also difficult to conduct. Therefore, molecular dynamic (M D) simulation is a useful tool to investigate the interfacial reinforcement mechanisms of graphene Œpolymer composites. In this chapter , molecular dynamics has been used to study interfacial properties of two nanocomposites, Graphene/polyethylene (PE) and C NT/polyethylene. Representative volume element (RVE) consists of a large numbers of amorphous polyethylene monomers that were modeled using Dreiding potential in LAMMPS (Hossain, Tschopp et al. 2010 ), while CHARMM potential were used for carbon -carbon interaction in CNT and Graphene (Stuart, Tutein et al. 2000). Periodic boundary conditions in three directions enabled us to simulate fiber pullout of a large graphene sheet. I n addition to pulling tests, tensile tests have been also performed on both 84 CNT/PE and Graphene/PE and elastic constants of these two composite systems have been compared. 4.2 Theory Molecular dynamics use s Newton™s equations of motion to calculate the trajectory of atoms in any system. The advantage of MD simulation is its capability of updating a system in future time step based on the current state of the system. At each step, the forces on atoms are comp uted and combined with current position as velocity to update positions and velocities of atoms in t smaller time steps by using these equations: = (1) = (,,,–) (2) Where , : Mass of atom : Position of atom : Force on atom (,,,–): Potential energy in N -atoms system 4.3 United -atom model United -atom approach has been used to represent the polyethylene where the monomeric units ŸC := CH2 of a olyethylene c hain within the matrix are treated as single spheres (Hossain, Tschopp et al. 2010 ). Bond relations are static here, so no breaking or forming of bonds within the polyethylene matrix is possible. 85 The Dreiding potential (Hossain, Tschopp et al. 2010 ) is used in this study . The Dreiding potential has four contributing terms; bond stretching, changes in bond angle, changes in dihedral rotation, a nd van der Waals nonbonded interactions. The force field and its respective parameters for the PE system are given in the Table 4.1. The total force field energy can be expressed as E = E bond + Eangle + Edihedral + Evan der waals (3) Where Ebond (r) = ½ Kb(r Œ r0)2 (4) Kb : force constant, r0 : equilibrium bond distance, r: bond distance Eangle () = ½ K (- 0)2 (5) K 0 Edihedral =30(cos) iiic (6) Kd : force constant, n: periodicity, d =4 (7) This is a Lennard -Jones potential model. parameter, r: separation of two atoms For the Lennard -Jones potential equations, there is an additional switching function that ramps the energy smoothly to zero between an inner and outer cut -off therefore preventing discontinuities in the energy and force. 86 Table 4.1 . Parameters for the bond stretching, the bond angle bending, and the torsional or dihedral potentials to model the intramolecular interaction of a linear polyethylene chain in the framework of an united -atom approach (Hossain, Tschopp et al. 2010 ). Interaction Form Parameters Bond Length Ebond (r) = ½ Kb(r Œ r0)2 Kb=350 kcal/mol , r 0=1.53 A Bond Angle Eangle () = ½ K (- 0)2 K=60 Kcal/mol/rad 2 , 0=1.911 rad Dihedral Angle Edihedral =30(cos) iiic C0=1.736, C 1=-4.49, C2=0.776, C3=6.99 (Kcal/mol) =4 CHRAM potential has been used to model CNT, Graphene sheet and xGnP and also the interaction between these nanoparticles with PE has been modeled using this potential (MacKerell, Bashford et al. 1998 , Stuart, Tutein et al. 2000 ). Dreiding and CHARM potential have been summarized i n Figure 4.1. 87 Figure 4.1 . Summary of potential s used for CNT, Graphene and polyeth ylen e 4.4 Fundamental algorithms in MD simulation MD simulation primarily concerns about solving Newton™s equations of motions for N - atoms system. At each timestep the total energy or Hamiltonian of the system is written as: 111 111(,,,...)(,,...) NN HKppppErrrr (8) Etotal Ebonded Enonbonded Ebond Eangle Edihedral Eelectrostatic EvanderWaals Dreiding Dreiding CHARMM 20() EK ikelectrostatic nonbonded ik pairs qqEDr126 ikik LennardJones nonbonded ikik pairs ACErr0(cos) iiiEC1,4 (1cos()) pairs EKn 201,2 ()bondstretch bpairs EKbb 20()UBUBEKSS (Dreiding) (CHARMM) 88 21111(,,,...) 2NiNiipKpppp m (9) ip : Momentum of i-th atom Based on the total energy conservation criteria the equations of motion of the N -atoms system can be defined in terms of a system of ordinary differential equations: 0iiiprm (10) 12,3,... (,) 0NiiErrrr pr (11) Let™s assume the above system of equations is obtained at time t. Finite difference method is used to update the current trajectory after solving the equations of motion. At the current state, acceleration of each atom ()ii atr is calculated by combining the Equations ( 10) and (11) as: 12,3,... (,) 0NiiiErrrr mrr (12) This acceleration is combined with positions and velocities at time t to calculate the positions ()irt and ()ivt velocities of the atoms at time t tand so on. This is commonly known as fiTime integration (TI) schemefl in MD simulation. There are several types of TI schemes available. Verlet algorithm is the most widely used among all the available methods. This algorithm uses the positions and accelerations at time t - t to predict the position and velocities at time t t. The advantage of Verlet algorithm is its simpl icity , speed and 89 memory efficien cy. So, Verlet algorithm i s highly used in large systems. The positions are calculated at t t by using velocities and accelerations at t by using the following equation : 21()()()() 2iiii rttrttvttat (13) The velocities at time 1()2tt is calculated using the velocities and accelerations at is calculated using the velocities and accelerations at time t: 11()()() 22iivttvttat (14) Hence, using these velocities and positions the forces on each atom is computed at current positions. This gives ()iatt . Finally, the velocities at tt are calculated: 151 ()()()()() 366 iiiii vttvttattattatt (15) 4.5 Ensemble average in MD simulation As discussed earlier the main objective of MD simulation is to replicate a material™s structure at atomic scale and predict macroscopic properties without directly depending on experiments. A large system of atoms (approximately 1023 ) might represent a system equivalent to macroscopic level. At each timestep t the force as well as acceleration on each atom will be calculated. Finally, any property average of this large system can be calculated by averaging for very long time t : 0lim[(),()] Average ii tptrtdt (16) [(),()] ii ptrt : Function of if (17) 90 There are different types of ensemble averages used in MD simulation. Among these first three types of ensembles are widely used in MD simulations. Based on different conditio ns these are: - Constant energy, constant volume (NVE) - Constant temperature, constant volume (NVT) - Constant temperature, constant pressure (NPT) - Constant temperature, constant stress (NST) - Constant pressure, constant enthalpy (NPH) In the current work the systems were equilibrated by using both NPT and NVT conditions. 4.6 Nanoscale RVE A single -walled CNT , Graphene sheet or xGnP embedded into an amorphous polyethylene matrix comprised the unit cell s of the MD simulation s. The CNT/PE , Gr/PE or xGnP/PE composites were created according to the following procedure. Initially, the PE amorphous structure has been create d; the coordinates of chains united atoms were randomly generated on an FCC lattice, using a random walk method similar to that of Shepherd (Shepherd 2006 , see Figure 4.2) . Then we o btain ed the equilibrium configuration of amorphous PE; all the chains are equilibrated for 106 MD time steps in the isothermal -isobaric (NPT) ensemble at 300 K to allow for energy relaxation and sufficient entanglement. In the case of CNT a cylindrical pore should be generated ; a cylinder indenter is placed in the center of the simulation box which exerts a forc e of magnitude F(r) = k(r R)2 on each united atom, where k is the specified force constant, r is the distance between the united atoms and the cylinder axis x, R =R CNT+Rcutoff , RCNT is the 91 radius of the CNT, and Rcutoff = 2.5 is the van der Waals cutoff distance. Finally a CNT has been embeded into the cylinder pore; equilibrate the CNT/PE composite for 500,000 MD time steps at 300 K using NPT dynamics to create a zero initial stress state. The same procedure has been repeated for Graphene/PE and x GnP/PE systems and RVEs are shown in Figure 4.3. The CNT , Gr and xGnP spanned the total length of the unit cell s, as shown in Figure 4.3. Periodic boundaries were applied in x and y directions. The t wo edges of the unit cell in the z direction were fixe d to prevent the rigid body motion of the PE . Figure 4.2 . Generating PE chains using a self -avoiding walk algorithm Avoiding Walk Density=0.1 united atom/volume Density=0.3 united atom/volume 92 Figure 4.3. (a) Nano RVE for CNT/PE, (b) Nano RVE for Gr/PE and (c) nano RVE for xGnP/PE 21.3 nm 20.6 nm 21.1 nm 9.1 nm CNT Graphene Exfoliated Graphite Nanoplatelets (xGnP), stacked Graphene a b c 93 4.7 Measuring interphase thickness for Graphene/PE system Fig ure 4.4 is extracted from the snapshot of the final MD simulated structure s, showing the carbon nanotube and graphene platelet polyethylene with interphase . One of the main characteristics of the interphase is its thickness , which can be defined as the distance from which the atomic properties (atomic density profile, atomic mobility) is different from those in bulk systems. The atomic concentration profiles of polymer (Fig ure 4.5) show that the polymer matrix reach its uniform bulk density beyond 3 5 Angstroms of graphene surface. Thus, we can estimate the interphase thickness as 3 5 nm. As part of the future work, graphene nanoplatelet will be modeled surroun ded by an interphase layer. The mechanical properties of this thin layer of polymer around the nanoplatelet also will be measured using molecular dynamic simulation developed in this chapter. Figure 4.4 . Polymer chains around the fillers in equilibrium state 94 Figure 4.5. Concentration profiles of polymer around Graphene sheet 4.8 Comparing interfacial properties of CNT and Graphene in polymer In this section molecular dynamics simulations of carbon nanotube (CNT) and Graphene pull -out from a polymer matrix were carried out. As the CNT or Graphene pull -out develops, variations in the displacement and velocities of the CNT and Graphene are monitored. The existen ce of a carbon ring -based period in Graphene and CNT sliding during pull -out is identified. Linear trends in the CNT and Graphene velocity Œforce relation are observed. The entire CNT and graphene pull -out process is simulated by MD. As the sample cell is periodic, the pull -out is simul ated by pulling the nanoparticle through, rather than out of the PE matrix. A unidirectional force is applied to each atom of the nanoparticles along the axial axis. The applied force is increased incrementally over time (Fig ure 4 .7). The simulation is run for ~100,000 MD time steps of 0.5 fs. The MD simulation captures the entire pull -out process. 00.5 11.5 22.5 050100 Relative Concentration Distance from Graphene Surface Polymer 95 Initially, sliding begins after a sufficient amount of force is applied to the nanoparticle . In Fig ure 4.7, the increments of appl ied force are shown, and over time, the velocity and the displacement (Figure 4.8) of the nanoparticles are monitored. In the simulation, the axial velocity of the nanoparticles is sufficiently greater than the velocity of the surrounding polymer, and the latter is therefore neglected in the analysis. Figure 4.6 . Displacement of individual carbon atoms of Graphene during the pull -out Timesteps(0.5fs)GrapheneDisplacement(A)050000100000150000200000-2-101234567891011122.5 A 96 Figure 4.7 . Increments of axial force applied to the nanotube and Graph ene over time. Velocity of the center of mass of a carbon NT and Graphene during the simulated pullout . Figure 4.8 . Displacement of individual carbon atoms of Graphene and CNT during the pull -out process Timestep(0.5fs)Velocity(A/fs)Totalforce(N)0100000200000300000-0.2-0.100.10.20.30.40.50.60.70.801E-072E-073E-074E-07ForceGrapheneVelocityCNTVelocityx10-3Timestep(0.5fs)Displacment(A)0100000200000300000051015202530CNTGraphene97 MD simulations of uniaxial tensile test on the Nanoscale RVEs (Gr/PE) shown in Figure 4.3 also have been performed. As it is shown in Figure 4.9, nanocomposite with 10 layers of graphene has higher modulus. Figure 4.9. Axial tensile responses of Graphene+PE (1,5 and 10 Graphene layers) 4.9 Conclusion For identical polymer systems (number of chains, number of monomers per chain) and strain rate, results show that the interfacial strength is lower in Graphene/PE system compared to CNT/PE. A chemically functionalized Graphene should be used to improve the interfacial bonding between the Graphene and the polymer. We expect that the chemical modifications of Graphene will increase the interfacial bonding characteristics between the Graphene and the polymer matrix. Hence, by modifying Graphene surface with ch emical surfactant, Graphene polymer composite can compete with CNT/polymer composite in interfacial and elastic 025507500.10.20.30.41200-50-101200-50-51200-50-1StrainStress (GPa)98 properties, while it maintains its superior electrical and thermal properties. As a future work t he effects of representative volume element siz e and polymer chain length and strain rate on these results also will be investigated. 99 REFERENCES 100 REFERENCES Hossain, D., et al. (2010). "Molecular dynamics simulations of deformation mechanisms of amorphous polyethylene." Polymer 51(25): 6071 -6083. MacKerell, A. D., et al. (1998). "All -atom empirical potential for molecular modeling and dynamics studies of pro teins." The Journal of Physical Chemistry B 102(18): 3586 -3616. Shepherd, J. E. (2006). "Multiscale modeling of the deformation of semi -crystalline polymers." Stuart, S. J., et al. (2000). "A reactive potential for hydrocarbons with intermolecular inte ractions." The Journal of chemical physics 112(14): 6472 -6486. 101 Chapter 5 Damage model for thermoplastic polymer nanocomposite In Chapter 3 elastoplastic material model has been used to model the deformation of HDPE and there was no damage in HDPE. In reality , both matrix cracking and fiber debonding occur during the deformation in nanocompoites. Here in this chapter , the damage progression in the matrix will be studied. 5.1 Damage model for high density polyethylene K46 -06-185 is a natural high molecular weight, high -density copolymer polyethylene developed for automotiv e fuel tanks. It combines excellent processability and outstanding physical performance, particularly environmental stress crack resistance (ESCR) and impact properties. This grade™s high melt strength aids the continuous extrusion production of multi - layer tanks. Properties of high -density copolymer polyethylene studied in this work have been listed in Table 5.1. High density polyethylene is a semi -crystalline thermoplastic polymer that can undergo large ductile deformation before losing all loa d-carrying capacity. The loss of load -carrying capacity results in the progressive degradation of material stiffness. The modified Gurson -Tvergaard -102 Needleman (GTN) model was applied to HDPE to model its plasticity and damage behavior. The GTN model was ori ginally developed to study damage evolution in porous media (Gurson 1977 , Tvergaard and Needleman 1984 ). Recently, t he original GTN model and its modified version s have been used to model the deformation of polymer s. For example, Oral et al. 2010 calibrated the original GTN model for polymer , and the result of their calibration suggested a different range of GTN parameters from that of a metal . Boisot et al. 2011 studied damage progression in nylon 11 using modified GTN model. Zaïri et al. 2005 studied the mechanical response of rubber toughened PMMA using original GTN model. Here in this study modified GTN model developed by Xue 2008 . has been used and calibrated to capture the damage progression in HDPE. Table 5.1 . Properties of K46 -06-185 (High density polyethylene copolymer) www.ineos -op.com Parameter SI Units ASTM Method Density 0.946 g/cc D4883 Melt Index 190°C/ 21.6 kg 4.2 g/10 min D1238 Compression Molded Samples D638 Tensile Strength (2 in/min) @ Yield 23.4 MPa @ Break 37.9 MPa Elongation (2 in/min) D638 @ Yield 12% @ Break >1000 % Flexural Modulus D790A Tangent Method 1040 MPa 2% Secant Method 807 MPa Notched Izod Impact Strength 68 kJ/m2 D256 Modulus of Elasticity 1035 MPa Poisson™s ratio 0.4-0.45 103 5.1.1 Modified Gurson -Tvergaard -Needleman (GTN) model The Gurson ŒTvergaard ŒNeedleman (GTN) porous plasticity model is a coupled damage criterion which is based on the effects of micro -void growth. In coupled ductile damage criteri on, damage accumulation is integrated in the constitutive equations by the evolution of the yield surface of the materials according to changes in some damage -induced density. In the GTN model, the behavior of the void ed solid is defined by pressure -sensitive plastic flow. The void volume fra ction is considered as the damage variable in the constitutive equation. In the GTN (Gurson -Tvergaard -Needleman) model , the flow potential is defined as, (1) Where is the yield function, is the hydrostatic pressure, is the prescribed hardening law for the matrix material as a function of the accumulated effective strain, , while , and are material parameters. The flow potential in Eq. ( 1) can be rewritten as, . (2) The effective porosity, , is determined by the following relationship: , (3) 2*2213()3(,,) 2cosh10 2YMMM fqpffq qf ()Yf11()33kkptr ()MMM ()MMd1q2q3q*22313(,,)()12cosh 02MYMMqpffqffq f*1,1,cccccfcfff fffff ff qff 104 Where the parameter represents the porosity, the constant is the porosity to trigger coalescence and the parameter represents the porosity at fracture. Note that the effective porosity after coalescence activation in Eq. (3) is differently defined in the ABAQUS embedded material as: , (4) Where . (5) In the GTN model, the evolution of the porosity is given by the sum of both nucleation and growth mechanisms, as: , (6) Void nucleation occurs by the cracking and decohesion of second -phase particles in ductile material s. The void nucleation is not homogenous due to the complexity of the microstructure. Chu and Needleman suggested a plastic strain controlled nucleation in which (Chu and Needleman 1980 ) . (7) Here, , and are the volume fraction of the nucleated voids, mean strain for void nucleation and its standard deviation, respectively. Void nucleation rate NMdfd as a function of the equivalent plastic strain of the undamaged matrix has been plott ed in Figure 5.1. The values of other parameters such as f N, SN, N have been obtained through the calibration procedure explained in section 5.1.2. fcfff*,(), cfcccc fcfff ffffffff ff 2113 3fqqq fq NGdfdfdf 21exp 22NNNMNMMNNfdfdAd SS NfNNS105 Figure 5.1 . Void nucleation rate The matrix is assumed to be incompressible, but the macroscopic response of the material is compressible due to the voids. So, df G is caused by the total volume change. The void growth mechanism is defined by, (8) Where is the volume increment of the plastic strain, . Figure 5.2 shows the contribution of void growth and void nucleation in porosity as a function of the plastic strain in the uniaxial tension test. As it was expected void nucleation reaches to the constan t value around strain equal to 0.2, consistent with Figure 5.1 where the Void nucleation rate NMdfd=0. Void growth keeps increasing until failure. Equivalent plastic strain of the undamaged matrix M0.0 0.1 0.2 0.3 0.4 dfN/d M0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 Void nucleation rate eN = 0.11 (1):(1) Gpp kkdffdfd Ipkkdpd106 Figure 5.2 . Void evolution during uniaxial tension Figure 5.3 also shows the contribution of void growth and void nucleation in porosity as a function of the plastic strain in the pure shear test. Void growth is zero during the shear deformation due to the lack of hydrostatic pressure in the shear mode. I t means if the total volume fraction to be nucleated is less than the void volume fraction at fracture, then the material is predicted not to fail at all . This is one of the shortcoming of the original GTN model. Xue et al. developed a modified GTN model t o overcome this shortcoming. Here in this study the ir modified GTN model has been used to model the damage in HDPE. Plastic strain 11 0.0 0.1 0.2 0.3 0.4 Porosity f 0.00 0.02 0.04 0.06 0.08 0.10 0.12 Void nucleation Void growth Todal voids 107 Figure 5.3 . Void evolution during pure shear Xue™s shear mechanism (Xue 2008 ) The s hear mechanism proposed by Xue 2008 was derived from geometrical relations of a cell structure with a circular void at the center subjected to a simple shear strain . During the shear deformation an initially spherical void become s spheroidal and the free surfaces of voids move closer to the boundary of the void containing cell . This is called void shearing. This has been shown in Figure 5.4. Since the material is incompressible under shear due to lack of the hydrostatic pressure ch ange this portion of the damage is missing in GTN model. So an artificial strain has been introduced by Xue 2008 to take into the account the damage of void shearing. 108 Figure 5.4 . Void shearing (Weck, Wilkinson et al. 2006 ) A 2D void shearing mechanism has been shown in Figure 5. 5. In this problem, a sq uare cell with length L contains a cylindrical void of radius R at the center, under the simple shear straining, the void elongates in the preferred direction. The relative position of the void does not change with respect to the cell due to the volume con servation of the cell . However, t he minimum distance between the free surface of the void and the boundary of the representative volume element decreases as the shear strain increases. Figure 5.5 . A schematic drawing illustrates the void shear mechanism. (Xue 2008 ) 109 The initial minimum distance between the free surface of the void and the cell boundary is 2LaR (9) In a simple shear strain, we have tan (10) Where is the deformation angle shown in Fig . 5.5. The minimum distance at a deformed configuration becomes 21cos 1aaa (11) It should be noted that this minimum distance is changing with the plastic flow. Using the logarithm definition of strain, an fiartificialfl strain associated with the reduction of the minimum distance is define d, 2loglog1 rot aa (12) For small , this strain reduces to 212rot (13) Comparing this artificial strain with the maximum value at which the free surface of the void reaches the boundary of the cell, the relative damage associated with the shearing of the void is 2log1 log 2rot DLR (14) Using Taylor™s expansion, the above Eq. can be approxima ted, i.e. 110 21214rot Df (15) Knowing =3 eq for simple shear and for small void volume fractions, the above Eq. can be approximately written as (1/2)2 3rot eqDf (16) For the 3D case, similar relation can be obtained, i.e. (1/3) (1/3)2 362rot eqDf (17) The incremental void shear damage will be 43qrot eqeq dDqfd (18) Where q3= 3.39 and q4 = 1/2 for 2D or q3 = 3.72 and q4 = 1/ 3 for 3D problems respectively for the present idealized cell structure. In the present study, Due to the complexity of the actual material structures, these constants are material dependent and have been calibrated by fitting to the experimental results. But q3 = 3.72 and q4 = 1/3 have been used as initial values for q3 and q4 in our calibration. Lode angle dependent GTN model Equation ( 18) is for the simple shear situation. In general loading case s, the damage relationship becomes complicated. For the shear mechanism in the modified GTN model, . (19) Here, and are material parameters, while is the Lode angle function defined as, 540qSMMdfqfgd 4q5q0g111 , (20) Where the normalized Lode angle, , is defined as, (21) Here, , and are principal deviatoric stress components. Under the plane stress condition, for the simple tension condition, for the pure shear condition, for the balanced biaxial tension condition. The shear mechanism described above can be easily added to the GTN model. This is undertaken by introducing each mechanism on the rate of the void volume fraction, f. Thus, the evolution of the damage variable can be rewritten as: NGS dfdfdfdf (22) The plastic strain increment in the GTN model is, , (23) Which is dependent on hydrostatic pressure, ; i.e., , (24) Where is the deviatoric stress. From the plastic work equivalence principle, (25) Therefore, 01g 21313266arctan ,11 3SSS SS 1S2S3S1 01pddp212*223133sinh 22312cosh 2MYMqpfqq fpppqpqffq SS1:3p SII:(1) pppYMMdwdfdfd 112 (26) Or (27) The linear elastic constitutive law for the Jaumann (or co -rotational) stress rate of the Cauchy stress is, . (28) Here, , and , are the elastic stiffness tensor, the total strain increment and the elastic strain increment, respectively. Then, the potential criterion described in Eq. ( 1) leads to, (29) From Eq. ( 27), Eq. ( 30) becomes, , (30) Where . (31) :(1) MMddf(1) :MMfddepdddddd CCC Cded5***40*(,,): :(1): 0MMMMMMMMMqNMMMfdfd ddfff ddd fAdfdqfgd ff CI**::(1) 0MMMM ffdddd fffd C540::(1): (1) (1) qNMMMMMddfdf Afqfg dddf f I113 After further manipulations from Eq. ( 31), becomes, (32) Where , (33) (34) And from Eqs. (3) and (4), (35) or . (36) Substituting Eq. ( 32) into Eq. ( 28), . (37) d**:::(1) MMMM ddfdf fffd CC*2231312cosh 2MMqpqffq *231**22313cosh 2312cosh 2MMMqpqfq fqpqffq *11,11 ,cccfcff ffff dfqff *1,,cfcccfcff ffffff fff**:: ::(1) MMMM fdf fffd CCCC114 5.1.2 Calibration of modified GTN model for high density polyethylene In this section material parameters for the modified GTN model will be calibrated for HDPE. Application of GTN model for polymers is limited mainly because the void nucleation, growth, and coalescence process is not present in most polymers in the same way it is in metals. However, the GTN model is a continuum softening and failure model, so by choosing an appr opiate set of parameters , this model will be suitable for modeling the failure of polymers (Imanaka, Nakamura et al. 2003 , Zaïri, Aour et al. 2008 , Zaïri, Naït -Abdelaziz et al. 2008 , Oral, Anlas et al. 2010 ). True stress -strain curves extracted from uniaxial tensile test were utilized for the calibration of elastic properties (see Figure 5.6). Young™s modulus E=1018.55 MPa a nd Poisson™s ratio =0.4 0, which is in a good agreement with the information in Table 5.1. Figure 5.6 . Stress -strain response (uniaxial tensile test) head speed=0.8 mm/min 0510 15 202530 35 40 4500.1 0.2 0.3 0.4 0.5 True stress (MPa) True strain HDPE 115 Plastic strain histories were measured by eliminating elastic strains from total strain obtained from DIC (digital image correlation) tests. Measured plastic strain histories along 22 d irection (transverse direction) are dependent on the assumed Poisson™s ratio (Figure 5.7) . Measured plastic strain ratio deviate from the isotropic curves as plastic strain 11 increases, which implies that void nucleation and growth has occur red. Figure 5.7. Plastic strain histories as a function of Poisson™s ratio compared with undamaged material (isotropic material) The m odified GTN model has been calibrated for HDPE in the following steps: (1) From the measured plastic strains during the uniaxial tension test Ł 6 GTN parameters ( q1, q2, SN, fN, eN and f 0) for the original GTN model excluding hardening parameters has been calibrated by utilizing the analytical code (Matlab) (2) From the measu red stress -strain curves Plastic strain 11 0.00 0.02 0.04 0.06 0.08 0.10 0.12 0.14 0.16 0.18 Plastic strain 22 -0.08 -0.06 -0.04 -0.02 0.00 Exp. (nu=0.45) Exp. (nu=0.44) Exp. (nu=0.40) Isotropic material (von Mises) 116 Ł 5 parameters for the hardening evolution of the undamaged HDPE can be calibrated by utilizing the analytical code (Matlab) (3) From the measured plastic strains during the simple shear test (using the notched specimen) Ł 2 additiona l parameters (q4 and q5) for the modified G TN model can be calibrated by performing FE simulations (ABAQUS and LS -OPT softwares ) Figure 5.8 shows the comparison of the calculated and measured plastic strain history. The original GTN model parameters have been obtained using this information as it is shown in this Figure. Figure 5.8 . GTN model, isotropic material (Von Mises) and experimental uniaxial tests Plastic strain 11 0.00 0.02 0.04 0.06 0.08 0.10 0.12 0.14 0.16 0.18 Plastic strain 22 -0.08 -0.06 -0.04 -0.02 0.00 Exp. (nu=0.40) G-T-N model Isotropic material (von Mises) 21231 01.2, 0.8, 0.06, 0.12, 0.11 0.0005NNN qqqq Sfe f 117 Limit of the porosity ff at fracture From equation 2, l imit of the porosity ff at fracture can be calculated. (38) The maximum value of ff is (39) Which is a function of q1 , q2, q3 and stress triaxiality ()Mp where The maximum allowable value for porosity (ff) at fracture as a function of stress triaxiality is plotted in Figure 5.9. The maximum allowable value for porosity (ff ) at fracture is maximum in pure shear and minimum in the case of balance biaxial mode. For different set of q1 and q2 , the shape of the curve is the same, but it shifts up or down based on the values of q1 and q2. For the balanced biaxial mode *axMfis 0.55. *2231312cosh0 2Mqpqffq 22213*ax333coshcosh 22MMMqpqpqqfq : Hydrostatic pressure : von Mises effective stress Mp212 23*ax32123 coshcosh 1, 0.55176 (when 1.2, 0.8, ) MMlqqqq pfqqqqq 118 Figure 5.9 . Limit of the porosity ff at fracture as a function of stress triaxiality Figure 5.10 shows the maximum allowable value for porosity (ff) at fracture as a function of q1 in the balanced biaxial mode . For small q1 (less than 1) the allowable porosity will drastically increase and as q1 becomes larger (more than 3) this value become s very small. Figure 5.11 shows the maximum allowable value for porosity (ff) at fracture as a function of q2 in balanced biaxial mode . Stress triaxiality -1.0 -0.5 0.0 0.5 1.0 Maximum porosity f* 0.0 0.2 0.4 0.6 0.8 1.0 q1=1.2, q 2=0.8, q 3=q 12119 q10.0 0.3 0.6 0.9 1.2 1.5 1.8 2.1 2.4 2.7 3.0 Maximum porosity f* 0.0 0.2 0.4 0.6 0.8 1.0 q2=0.6 q2=0.7 q2=0.8 q2=1.0 q2=1.2 Figure 5.10 . Limit of the porosity f f at fracture for the balanced biaxial mode , effect of q1 q1123Maximum porosity f* 0.0 0.1 0.2 0.3 0.4 0.5 0.6 120 q20.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 Maximum porosity f* 0.0 0.2 0.4 0.6 0.8 1.0 q1=0.6 q1=0.9 q1=1.2 q1=1.5 q1=2.0 Figure 5.11 . Limit of the porosity ff at fracture for the balanced biaxial mode , effect of q1 121 Figure 5.12 . uniaxial tensile tests, GTN model versus experimental results The second step is obtaining the hardening parameters for the undamaged matrix. Curve fitting based on combined Swift -Voce hardening law was carried out (Swift 1952 ). (40) Void coalescence trigger and void at fracture parameters, fc and ff , were assumed as 0.32 and 0.4, respectively. The results of the calibration of GTN model for HDPE have been summarised in Table 5.2 and 5.3. Engineering strain 11 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 Engineering stress 11 (MPa) 051015 20 2530Exp. G-T-N model (nu=0.4) Uniform deformation limit c48.535 0.98999 0()(1)72.102MPa(0.055637)19.860MPa(1) MMnMMMKBe e 122 Table 5.2 . Hardening evolution of the undamaged material K 0 n B c 72.102MPa 0.05564 0.989 19.8MPa 48.53 Table 5.3 . Material coefficients for modified G -T-N model q1 q2 q3 q4 q5 SN fN eN fc ff finit 1.2 0.8 1.44 2.65 0.42 0.06 0.12 0.11 0.32 0.4 .0005 The third step is calibration of modified GTN model from the measured plastic strains during the simple shear test (using the notched specimen in section 2.7). The parameters, q4 and q5 have been calculated by performing FE simulations with Abaqus and LS -OPT software. A python scripts has been written to automatically extract the load -displacement from FE simulation and compare it with experiment al data until the curves were matched. Figure 5.13 shows a n FE model of the shear specimen designed in section 2. 7. A very fine mesh has been used around the notche areas to guarantee the numerical stability of the modified GTN model. A total of 19782 hexahedral (C3D8R) elements has been used in this simulation. The boundary conditions are applied to restrict all fr eedom of movement at one end of the model and the other end is coupled to a reference point, which is further given a displacement outward; causing elongation. 123 Figure 5.13 . Mesh of shear specimen The result of calibration using Abaqus and Ls -opt software is shown in Figure 5.14. Damage evolution or void evolution during the deformation is shown in F igure 5.15. Initial void density was 0.0005 as obtained through calibration shown in the first frame of the figure. Void evolves within the notche area as it was expected. Figure 5.1 4. Load -displacement curve for shear specimen,experiment versus simulation 04080120 160 200 0246810Load (N) Displacement (mm) shear test, simulation Shear test (expriment) 124 Figure 5.15 . Void e volution during loading in the modified specimen 125 5.2 Damage model for HDPE/GNP nanocomposite 5.2.1 Representative volume element In this study, for each RVE model, at least 40 GNP were randomly distributed in the cubic matrix following the procedure explained in section 3.2.1. Since the GNPs are randomly distributed in t he matrix, the nanocomposites behave as an isotropic solid at the macroscale. The damage and fracture behavior are analyzed through such RVE models. A cubic RVE model and the finite element mesh of a polymer/ GNP nanocomposites generated in Abaqus software is shown in Figure 5.16. Due to the complexity of the local geometry, the matrix and GNPs were meshed with fir st-order tetrahedral elements ( C3D4 elements ). The RVE model is subjected to quasi -static uniaxial tensile loading to study its damage behavior wi th special consideration for the microstructure of actual polymer/ GNP nanocomposites. Figure 5.16 . RVE model with 1% weight fraction of GNP, Finite element meshes of (a) RVE (b) GNP a b 126 5.2.2 Boundary Conditions Bound ary cond itions are very crucial in analyzing an RVE to ensure the compatibility of deformation and precise calculation of stress and strain. Here in this study, simulation of the RVE under uniaxial stress was carried out and compared with actual experimental results . Displacement boundary conditions were applied to the surface of the RVE finite element mesh (see Figure 5.17) . (0,,)0 uyz (41) (,0,)0 vxz (42) (,,0)0 wxy (43) (,,) yy vxLz (44) Where, y is the prescribed displacement for uniaxial loading along the Y axis; x and zare computed displacements from the condition that the stress resultants on the surface ,0, xzxLxzL are equal to zero, as: 00(,,)0 yz LLxxLyzdydz (45) 00(0,,)0 yz LLxyzdydz (46) 00(,,)0 xy LLzzxyLdxdy (47) 127 00(,,0)0 xy LLzxydxdy (48) Where, (,,),(0,,) xxx Lyzyz and (,,),(,,0) zzz xyLxy are the stress distributions acting on the surface x=L x, x=0, z=L z and z=0. These zero resultant force conditions were enforced using multipoint constraint equations in Abaqus (linear constraint equations) (Abaqus 2012 ). The macroscopic nominal stress and strain can be obtained from: ,RyyxzFLL (49) yyyL (50) Where, ,RyFis the reactio n force on the top surface Figure 5.17 . Multipoint constraint boundary condition L L L 128 Explicit FEM has been used to solve the failure and damage problems in order to overcome calculation convergence issues. In the next section, the dynamic explicit formulation will be discussed. 5.2.3 Results The results of the simulation are shown in Figure 5.18. The damage always initi ates from the GNP platelets and then propagate inside the polymer. As can be seen, some parts of the polymer around the GNP have been partially damaged , but the rest of GNP can still effectively transfer load between the GNPs and the matrix. Figure 5.18: Damage evolution in GNP/HDPE composite. (SDV2: void density) 129 Homogeni zed engineering stress and strain of the RVE has been calculated using equations 49 and 50. Figure 5.19 compare s the engineering stress -strain obtained from simulation with the one from experiment. Figure 5.1 9. Overall engineering stress -strain for GNP/HDPE composite versus HDPE 5.3 Conclusion and future work A damage model of HDPE/GNP nanocomposite has been proposed to capture two damage mechanism s, matrix cracking and fiber debonding. A 3D microstructure model was incorporated in a damage modeling problem in a nanocomposite where damage initiation and damage progression have been modeled using cohesive -zone and modified Gurson -Tvergaard -Needleman (GTN) material models. To the best of author knowledge , this is the first 3D RVE mod el which takes into account two damage behavior s for nanocomposites 0510 1520253035020406080Eng Stress Eng. Strain HDPE+2.2 wt% GNP (exp) HDPE+2.2 wt% GNP (Simulation) HDPE (exp) 130 using modified GTN model and cohesive surface. This proposed hyb rid model, which was the goal of this project, can be used to model the damage initiation and progression of any type of nanoparticles and thermoplastic. It should be noted that there are limitation s to this model and further refinement are need ed to improve the accuracy of the simulation results. Further experimental validation and a through parametric study are necessary for such improvement , but are beyond the initial scope of this work. These are proposed as future work on this study: Modified tensile specimen to further calibration of modified GTN model and validation. The calibration procedure in this study should be validated using the modified dogbone specimen where there is a range of different stress triaxialities in the specimen. In thi s case the uniqueness of calibration parameters will be guaranteed. Study the effect of the interphase In this study, it was a ssumed that GNP is a single layer particle, however the interphase zone created by alter ed polymer chains in the vicinity of the nanoplatelets as it was simulated using molecular dynamics simulation in chapter 4 plays an important reinforcing role. The interphase has the same length -scale as that of the nanoplatelet and its properties were repo rted to be different (mostly higher) from those of the bulk polymer matrix (Shelley, Mather et al. 2001 , Sheng, Boyce et al. 2004 , Chen, Evans et al. 2008, Sikdar, Pradhan et al. 2008 ). The increase in stiffness was found to be due to changes in the mobility of polymer chains adjacent to nanoplatelet, as well as changes in crystallinity in semicrystalline polymers like Nylon -6 (Chen et al., 2008). A three phase composite comprises of GNP, polymer and interphase layer should 131 be constructed and its damage behavior should be modeled by utilizing the damage model developed in this study. Damage model for fiber reinforced polymer/GNP nanocomposites Polymer nanocomposites can not compete with continuous fiber reinforced polymers (FRPs) in strength and the modulus. A hybrid material consists of continuous fiber and polymer nanoparticle composite can be developed where the polymer/nanoparticle nanocompo site is treated as a nano -resin matrix (Qin, Vautard et al. 2015 ). Using the damage model for polymer nanocomposite presented here , a damage model will be developed to predict the initiation and progression of hybrid thermoplastic/ FRPs nanocomposites. Viscoplastic damage model for polymer Polymers are rate dependent. It means their properties vary by changing test speeds and also their properties are temperature dependent. The modified GTN model in this study should be extended by adding viscoplastic constitutive equations in this code to c apture the rate -dependence of the elastic, plastic and failure behavior of polymers in nanocomposite. Parameter study The p aramet ric study shoul d be performed to study the effect s of nanoparticle orientation, aspect ratio and volume fraction on the macroscopic mechanical properties similar to the one presented in chapter 3. 132 REFERENCES 133 REFERENCES Abaqus, A. (2012). "Standard: user™s manual version 6.12." Simulia, USA . Boisot, G., et al. (2011). "Experimental investigations and modeling of volume change induced by void growth in polyamide 11." International Journal of Solids and Structures 48(19): 2642 -2654. Chen, B., et al. (2008). "A critical appraisal of poly mer Œclay nanocomposites." Chemical Society Reviews 37(3): 568 -594. Chu, C. and A. Needleman (1980). "Void nucleation effects in biaxially stretched sheets." Journal of Engineering Materials and Technology(Transactions of the ASME) 102(3): 249-256. Gurson, A. L. (1977). "Continuum theory of ductile rupture by void nucleation and growth: Part IŠYield criteria and flow rules for porous ductile media." Journal of engineering materials and technology 99(1): 2 -15. Imanaka, M., et al. (2003). "Fracture t oughness of rubber -modified epoxy adhesives: effect of plastic deformability of the matrix phase." Composites Science and Technology 63(1): 41 -51. Oral, A., et al. (2010). "Determination of Gurson -Tvergaard -Needleman model parameters for failure of a pol ymeric material." International Journal of Damage Mechanics : 1056789510385261. Qin, W., et al. (2015). "Mechanical and electrical properties of carbon fiber composites with incorporation of graphene nanoplatelets at the fiber Œmatrix interphase." Composit es Part B: Engineering 69: 335 -341. Shelley, J., et al. (2001). "Reinforcement and environmental degradation of nylon -6/clay nanocomposites." Polymer 42(13): 5849 -5858. Sheng, N., et al. (2004). "Multiscale micromechanical modeling of polymer/clay nano composites and the effective clay particle." Polymer 45(2): 487 -506. Sikdar, D., et al. (2008). "Altered phase model for polymer clay nanocomposites." Langmuir 24(10): 5599 -5607. Swift, H. W. (1952). "Plastic instability under plane stress." Journal of the Mechanics and Physics of Solids 1(1): 1 -18. Tvergaard, V. and A. Needleman (1984). "Analysis of the cup -cone fracture in a round tensile bar." Acta metallurgica 32(1): 157 -169. 134 Weck, A., et al. (2006). "2D and 3D visualization of ductile fracture. " Advanced Engineering Materials 8(6): 469 -472. Xue, L. (2008). "Constitutive modeling of void shearing effect in ductile fracture of porous materials." Engineering Fracture Mechanics 75(11): 3343 -3366. Zaïri, F., et al. (2008). "Steady plastic flow of a polymer during equal channel angular extrusion process: Experiments and numerical modeling." Polymer Engineering & Science 48(5): 1015 -1021. Zaïri, F., et al. (2008). "Modelling of the elasto -viscoplastic damage behaviour of glassy polymers." International Journal of Plasticity 24(6): 945 -965. Zaïri, F., et al. (2005). "Constitutive equations for the viscoplastic -damage behaviour of a rubber -modified polymer." European Journal of Mechanics -A/Solids 24(1): 169 -182.