Multiscale Modelling of Nano-Clay Filled Shape Memory Polymer Foams A dissertation accepted by the Faculty of Energy-, Process- and Bio-Engineering of the University of Stuttgart in partial fulfillment of the requirements of the degree of Doctor of Engineering Sciences (Dr.-Ing.) By M. Sc. Marwan A. Salman born in Baghdad, Iraq Supervisor: Prof. Dr. rer. nat. Dr. h. c. Siegfried Schmauder Co-Examiner: Prof. Dr.-Ing. Holger Steeb Date of oral examination: 15. November 2019 Institute for Materials Testing, Materials Science and Strength of Materials University of Stuttgart 2020 Multiscale Modelling of Nano-Clay Filled Shape Memory Polymer Foams Von der Fakultät für Energie-, Verfahrens- und Biotechnik der Universität Stuttgart zur Erlangung der Würde eines Doktors der Ingenieurwissenschaften (Dr.-Ing.) genehmigte Abhandlung Vorgelegt von M. Sc. Marwan A. Salman aus Baghdad, Irak Hauptberichter: Prof. Dr. rer. nat. Dr. h. c. Siegfried Schmauder Mitberichter: Prof. Dr.-Ing. Holger Steeb Tag der mündlichen Prüfung: 15. November 2019 Institut für Materialprüfung, Werkstoffkunde und Festigkeitslehre der Universität Stuttgart 2020 Erklärung über die Eigenständigkeit der Dissertation Ich versichere, dass ich die vorliegende Arbeit mit dem Titel Multiscale Modelling of Nano-Clay Filled Shape Memory Polymers Foams selbständig verfasst und keine anderen als die angegebenen Quellen und Hilfsmittel benutzt habe; aus fremden Quellen entnommene Passagen und Gedanken sind als solche kenntlich gemacht. Declaration of Authorship I hereby certify that the dissertation entitled Multiscale Modelling of Nano-Clay Filled Shape Memory Polymers Foams is entirely my own work except where otherwise indicated. Passages and ideas from other sources have been clearly quoted. 30.01.2020, Stuttgart Marwan Salman Acknowledgments First and foremost, praises and thanks to God, the Almighty, for letting me through all the difficulties. I have experienced Your guidance day by day. You are the one who let me finish my degree. I will keep on trusting You for my future. Thank you, Lord. I would like to acknowledge my indebtedness and render my warmest thanks to my supervisor Prof. Dr. rer. nat. Dr. h. c. Siegfried Schmauder, you have been a tremendous mentor for me. I would like to thank you for encouraging my research and for allowing me to grow as a research scientist. Your advice on both research as well as on my career have been invaluable. I would also like to thank my co-supervisor, Prof. Dr.-Ing. Holger Steeb, for letting my defense be an enjoyable moment, and for your brilliant comments and suggestions, thanks to you. A special thanks to my family. Words can not express how grateful I am to my mother and father for all of the sacrifices that you’ve made on my behalf. Your prayer for me was what sustained me thus far. I would also like to thank my beloved wife, Asmaa. Thank you for supporting me for everything, and especially I can’t thank you enough for encouraging me throughout this experience. Furthermore, great thanks go for my son, Yousif, and my small angel, Sarah, who have been the light of my life for the last six years and who have given me the extra strength and motivation to get things done. I would like to say thanks to my friends and research colleagues, Dipl.-Ing. Wolfgang Verestek, Dipl.-Ing. Vinzenz Guski, Dr.-Ing. Peter Binkele, Dr. rer. nat. Galina Lasko, Dipl.-Math. Stefan Küster, Dipl.-Phys. Martin Hummel, mag. ing. mech. Marijo Mlikota, Dr. rer. nat. Alejandro Mora, M.Sc. Dennis-Michael Rapp and M.Sc. Saeid Sajadi, for providing a motivating working atmosphere. Many thanks for your constant encouragement. Finally, I would like to thank the German Academic Exchange Service DAAD (Deutscher Akademischer Austauschdienst), for providing the funding which allowed vi me to undertake this research. Contents Contents vii List of Symbols xi Abstract xvi Zusammenfassung xix 1 INTRODUCTION AND REVIEW OF LITERATURE 1 1.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Literature survey . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 1.3 Outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 2 THERMOMECHANICAL BEHAVIOUR OF SHAPE MEMORY POLYMER FOAMS 11 2.1 Material description and experimental observations . . . . . . . . . . 11 2.2 Multiscale modeling approaches . . . . . . . . . . . . . . . . . . . . . 14 2.3 Constitutive models . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 2.3.1 Rubbery phase . . . . . . . . . . . . . . . . . . . . . . . . . . 17 2.3.2 Glassy phase . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 2.3.3 Shape recovery effect . . . . . . . . . . . . . . . . . . . . . . . 21 3 MOLECULAR DYNAMICS SIMULATION OF NANO-CLAY FILLED SHAPE MEMORY POLYMERS 23 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 vii Contents viii 3.2 Material description . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 3.3 Molecular Dynamics (MD) Method . . . . . . . . . . . . . . . . . . . 28 3.3.1 Statistical ensembles . . . . . . . . . . . . . . . . . . . . . . . 30 3.3.1.1 NVE ensemble . . . . . . . . . . . . . . . . . . . . . 30 3.3.1.2 NPT ensemble . . . . . . . . . . . . . . . . . . . . . 31 3.4 Model description and calculations . . . . . . . . . . . . . . . . . . . 38 3.4.1 Coarse-grain parametrization and Force-Field . . . . . . . . . 38 3.4.2 Model preparation . . . . . . . . . . . . . . . . . . . . . . . . 41 3.4.3 Calculation of Elastic Constants . . . . . . . . . . . . . . . . . 47 3.4.4 Calculation of the interfacial strength parameters . . . . . . . 48 3.5 Results and discussion . . . . . . . . . . . . . . . . . . . . . . . . . . 50 3.5.1 Elastic Constants . . . . . . . . . . . . . . . . . . . . . . . . . 50 3.5.2 Transition temperature and thermal expansion coefficient . . 51 3.5.3 Fracture characterization of the epoxy/clay interface . . . . . 53 3.6 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 4 MICROSCOPIC MODELLING OF SHAPE MEMORY POLYMER NANO-COMPOSITES 63 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 4.2 Homogenization and multiscale modeling approaches . . . . . . . . . 66 4.2.1 Multiscale finite element method . . . . . . . . . . . . . . . . 66 4.2.2 Homogenization methods . . . . . . . . . . . . . . . . . . . . . 67 4.3 Analytical models for heterogeneous materials . . . . . . . . . . . . . 69 4.3.1 Mori-Tanaka (MT) model . . . . . . . . . . . . . . . . . . . . 69 4.3.2 Self Consistent model . . . . . . . . . . . . . . . . . . . . . . . 72 4.4 Model preparation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 4.4.1 Representative Volume Element (RVE) generation . . . . . . . 74 4.4.2 Representative Volume Element (RVE) size optimization . . . 75 4.4.3 Periodic Boundary Conditions . . . . . . . . . . . . . . . . . . 78 4.5 Determinations of effective parameters . . . . . . . . . . . . . . . . . 83 Contents ix 4.5.1 Effective elastic parameters . . . . . . . . . . . . . . . . . . . 83 4.5.2 Uniaxial finite deformation numerical test of the rubbery and glassy phases . . . . . . . . . . . . . . . . . . . . . . . . . . . 85 4.6 Results and discussions . . . . . . . . . . . . . . . . . . . . . . . . . . 87 4.6.1 Elastic constants . . . . . . . . . . . . . . . . . . . . . . . . . 87 4.6.2 Nonlinear hyperelastic numerical test for high temperature . . 89 4.6.3 Nonlinear viscoelastic numerical test for low temperature . . . 91 4.7 Summary and conclusions . . . . . . . . . . . . . . . . . . . . . . . . 93 5 MESOSCALE II MODELLING OF SHAPE MEMORY POLYMER FOAMS 97 5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97 5.2 Theoretical concepts of cellular materials . . . . . . . . . . . . . . . . 99 5.2.1 Fundamental foam concepts . . . . . . . . . . . . . . . . . . . 99 5.2.2 Effective stiffness analytical models . . . . . . . . . . . . . . . 100 5.2.3 Terminology of the polymeric open cells foams . . . . . . . . . 102 5.3 Computational homogenization scheme . . . . . . . . . . . . . . . . . 103 5.4 Material structure at mesoscale II . . . . . . . . . . . . . . . . . . . . 105 5.4.1 Liquid-state vs solid-state foaming process . . . . . . . . . . . 106 5.4.2 Epoxy shape memory polymer foam structure . . . . . . . . . 107 5.5 Model preparation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107 5.5.1 Geometric construction of 3D open-cell foam structure . . . . 107 5.5.1.1 Realistic 3D mesh reconstruction . . . . . . . . . . . 109 5.5.1.2 Repetition of tetrakaidecahedron (Kelvin) unit cells . 111 5.5.2 Boundary Conditions . . . . . . . . . . . . . . . . . . . . . . . 112 5.5.3 Contact installation . . . . . . . . . . . . . . . . . . . . . . . . 114 5.6 Results and discussions . . . . . . . . . . . . . . . . . . . . . . . . . 115 5.6.1 Effective elastic and lame constants . . . . . . . . . . . . . . . 115 5.6.2 A numerical large deformation compressive test at high tem- perature . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117 List of Symbols x 5.6.3 A numerical large deformation compressive test at low temper- ature . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120 5.7 Summary and conclusions . . . . . . . . . . . . . . . . . . . . . . . . 123 6 MACROSCALE MODELLING OF THE SHAPE MEMORY POLY- MER FOAMS 125 6.1 Numerical homogenization scheme . . . . . . . . . . . . . . . . . . . . 125 6.2 Thermomechanical parameter identification of glassy phase . . . . . . 126 6.3 Thermomechanical parameters of the rubbery phase . . . . . . . . . . 127 6.4 Numerical model implementation . . . . . . . . . . . . . . . . . . . . 128 6.5 Results and discussions . . . . . . . . . . . . . . . . . . . . . . . . . . 130 6.5.1 Finite strain numerical compression test at low temperature . 130 6.5.2 Finite strain numerical compression test at high temperature . 130 6.5.3 Shape recovery response of the SMPFs . . . . . . . . . . . . . 132 7 Conclusions and future work 135 7.1 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 135 7.2 Future work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 136 Bibliography 139 List of Figures 155 List of Tables 161 List of Symbols αζ Volumetric expansion coefficient F i Force acting on the particle i p Momenta I Fourth order identity tensor D Rate of stretching Dp Shape change rate of the relaxed configuration F Deformation gradient tensor F p Deformation gradient of the relaxed configuration I Second order identity tensor L Velocity gradient Lp Velocity gradient of the relaxed configuration N Normalized tensor aligned with the deviatoric driving stress state TN Network stress that captures the effects of orientation-induced strain hard- ening T ∗′ Deviatoric portion of the driving stress state T g Cauchy stress tensor of the glassy phase T r Cauchy stress tensor xi List of Symbols xii V e Hencky strain W Spin δfail Damage evolution δinit Damage initiation ∆G Zero stress level activation energy ϵ0 Vacuum dielectric permittivity ϵ Depth of the potential at the minimum r λchain Stretch on every chain in the polymer network t Tractions H Hamiltonian equations of motion L Lagrangian X 2 Chi-square criterion L Langevin function µg,λg Lamé constants ν Poisson’s Ratio ρ Density σint Internal stress σ Stress tensor σ Zero-crossing distance for the potential τ Effective equivalent shear stress θeq Equilibrium value of the angle among three bonded atoms List of Symbols xiii G̃ Metric tensor h̃ Cell transformation matrix s̃i Scaled coordinate vector of the particle i A Width of the phase transition zone c Heat capacity per unit mass C∗ Effective stiffness tensor Cr Initial hardening modulus E Young’s Modulus Eangle Energetics for angle bending among three bonded atoms Ebond Potential of the bond stretching between two bonded atoms Ecoul Coulombic Interactions between all pairs of atoms Ekin Kinetic energy Enonbond Non-bonded interaction between all pairs of atoms Etorsion Torsional energy among four bonded atoms Etot Total energy EV DW Attractive dispersion energy fg Volume fraction of the glassy phase fr Volume fraction of the rubbery phase G Shear Modulus H Concentration factor h Softening slope List of Symbols xiv k Boltzmann’s factor Kθ Constant of the energetics for angle bending among three bonded atoms Kn Interphase stiffness Kr Constant of the potential of bond stretching between two bonded atoms MP Inertia parameter MT Temperature parameter N Number of rigid molecular units between entanglements Nf Physical system’s degree of freedom P Thermodynamic pressure q Charge of the atom r Distance between two bonded atoms req Equilibrium bond distance s Athermal shear stress sss Preferred state T Absolute temperature t Time Tg Glass transition temperature Tr Reference temperature U (x) Potential that characterizes the interaction between the atoms U Potential energy v Velocity of the atoms List of Symbols xv W Internal virial work of the atoms γ̇p Plastic shear strain rate γ̇0 Pre-exponential factor Abstract xvi Abstract Shape memory materials (SMMs) are defined as materials that can recover their undeformed configuration after applying an external stimulus like light, tempera- ture, magnetic field, etc. These materials are involved in many applications i.e. aerospace, biomedical and automotive. Among all SMMs, shape memory polymers (SMPs) emerged as the most promising materials. SMPs have received the attention of scientists because of many interesting properties, namely, lightweight structure, biodegradable, and cheap to manufacture. Furthermore, they can be deformed up to 200% of their original configuration keeping the ability to restore their original shape. Recently, material scientists developed SMP nanocomposites by reinforcing them with nanoinclusions (nano-clays, carbon nanotubes, etc.). This step increased the degree of interest of SMPs and expanded the range of applications in industry and medicine. Moreover, producing shape memory polymer foam nanocomposites led to a revolution in the shape memory materials field. The resulting extra lightweight materials can show the shape recovery phenomena under finite deformation with higher mechanical properties. Nano-clay filled shape memory polymer foams have heterogeneous structures at cer- tain levels of observation. For instance, the structure of the material consists of a matrix material with randomly distributed voids of different sizes at the macroscale level. At a lower level of observation, nano-clay layers are randomly distributed through the polymer matrix. Modelling of such multifunctional materials is a big challenge according to the fact that the mechanics of the underlying atomistic nano- and continuum mechanical microstructures have considerable impacts on the material behavior at the macroscale level. Additionally, many thermomechanical phenomena arise from the lower scale levels. Within this work, a full multiscale modelling approach able to simulate the ther- momechanical behaviour of the nano-clay filled epoxy shape memory polymer foam under finite strain and at atomistic, micro-, meso- and macro-scale levels is presented. In addition, the nonlinear hyperelastic and viscoelastic behaviour of the material is Abstract xvii investigated as well as the shape recovery behaviour. Such an approach is applied by dividing the length scale of the problem into four different levels. The first, and smallest one is the atomistic scale level, in which a set of molecular dynamics models are proposed to determine the thermomechanical properties of the nanocomposite materials. Furthermore, the interfacial strength between the epoxy polymer matrix and the layers is estimated in all possible layer orientations. Moreover, the effects of nano-reinforcements on the thermal expansion coefficient as well as the transi- tion temperature of the epoxy shape memory polymer are estimated using the MD method. The second level is the mesoscale I. A combination of the Finite Element method in conjunction with a Representative Volume Element (RVE) concept and a multiscale numerical homogenization scheme are employed to investigate nano-scale strengthen- ing effects on the thermomechanical properties of the shape memory polymer. More- over, analytical models are presented and used to compute the overall effective elastic constants of the RVE. On the other hand, the cellular structure of the SMPF is studied at the mesoscale II level. Real and artificial 3D geometries are generated and used in the FE method. In addition, a computational homogenization scheme is used to include the 3D RVE that is generated at the lower scale level in the analysis. Moreover, a numerical finite strain compressive test is applied to both phases (rubbery and glassy) of the SMPFs. In the next steps, all the estimated parameters of the mesoscale II level are accu- mulated and applied directly into the macroscopic constitutive models of both the rubbery and glassy phases. The thermomechanical performance of the nano-clay filled shape memory polymer foam is examined at temperatures lower and higher than Tg. Furthermore, two different thermomechanical loading cycles are employed to investi- gate the shape recovery response of the SMPF nano-composites. Last but not least, it is shown that multiscale modeling of SMPF provides a deep understanding of the material behaviours especially those which are originating from micro- and nano-scale material heterogeneous structures. Moreover, the herein ap- plied Multiscale Materials Modelling (MMM) is considered as one of the most accu- Abstract xviii rate simulation methods which is able to model available as well as future complex structures of heterogeneous nanocomposite materials. Zusammenfassung xix Zusammenfassung Formgedächtnismaterialien (Shape Memory Materials (SMMs)) sind definiert als Materialien, die ihre Ausgangskonfiguration nach Anlegen einer äußeren Anregung wie Licht, Temperatur, Magnetfeld usw. wiederherstellen können. Diese Materialien werden in vielen Anwendungen eingesetzt, wie z.B. in der Luftfahrt, der Biomedizin oder der Automobilindustrie. Von allen SMMs erwiesen sich Formgedächtnispolymere (Shape Memory Polymers (SMPs)) als die vielversprechendsten Materialien. SMPs haben aufgrund vieler interessanter Eigenschaften Beachtung gefunden, nämlich ihre leichte Struktur, ihre biologische Abbaubarkeit und ihre günstige Herstellung. Darüber hinaus können sie bis zu 200% von ihrer ursprünglichen Konfiguration verformt werden, wobei die Fähigkeit zur Wiederherstellung ihrer Ausgangsform erhalten bleibt. In jüngster Zeit haben Materialwissenschaftler SMP-Nanokomposite entwickelt, indem sie sie mit Nano-Einschlüssen (Nano-Tone, Kohlenstoff-Nanoröhrchen, etc.) verstärkt haben. Dieser Schritt steigerte das Interesse an SMPs und erweiterte das Anwendungsspektrum in Industrie und Medizin. Darüber hinaus führte die Herstellung von Formgedächtnispolymerschaum (Shape Memory Polymer Foams (SMPFs))-Nanokompositen zu einer Revolution im Bereich der Formgedächtnis- materialien. Die daraus resultierenden extra leichten Materialien weisen höhere mechanische Eigenschaften auf und zeigen das Phänomen der Formwiederherstellung bei endlicher Verformung. Nano-Clay gefüllte Formgedächtnis-Polymerschäume weisen auf bestimmten Beobachtungsebenen heterogene Strukturen auf. So besteht die Struktur des Materials auf der Makroebene beispielsweise aus einem Matrixmaterial mit zufällig verteilten Poren unterschiedlicher Größe. Auf einer niedrigeren Beobachtungsebene sind die Nano-Clay-Schichten zufällig im Polymer verteilt. Die Modellierung solcher multifunktionaler Materialien ist eine große Herausforderung, da die Mechanik der zugrunde liegenden atomistischen nano- sowie der beteiligte kontinuumsmechanis- chen Mikrostrukturen erhebliche Auswirkungen auf das Materialverhalten auf der Zusammenfassung xx Makroebene hat. Darüber hinaus resultieren viele thermomechanische Phänomene aus den unteren Skalen. Im Rahmen dieser Arbeit wird ein vollständiger multiskaliger Modellierungsansatz vorgestellt, der in der Lage ist, das thermomechanische Verhalten eines mit Na-nöton gefüllten Epoxid-Formgedächtnispolymerschaums unter äußerer Belastung auf atom- istischer, mikro-, meso- und makroskaliger Ebene zu simulieren. Darüber hinaus wird das nichtlineare hyperelastische und viskoelastische Verhalten des SMPFs sowie das Rückstellverhalten des Schaums untersucht. Dieser Ansatz wird angewendet, indem die Längenskala des Problems in vier verschiedene Ebenen unterteilt wird. Die erste und kleinste ist die atomistische Ebene, in der eine Reihe von Molekular- dynamikmodellen vorgeschlagen werden, um die thermomechanischen Eigenschaften der Nanokomposite zu bestimmen. Darüber hinaus wird die Grenzflächenfestigkeit zwischen der Epoxid-Polymermatrix und den Schichten für alle möglichen Schichtori- entierungen abgeschätzt. Zusätzlich werden die Auswirkungen von Nanoarmierungen auf den thermischen Ausdehnungskoeffizienten sowie die Übergangstemperatur mit der MD-Methode abgeschätzt. Die zweite Ebene ist die Mesoskala I. Eine Kombination aus Finite-Elemente- Methode, Repräsentatives-Volumen-Element(RVE)-Konzept und numerischem Homogenisierungs-schema wird eingesetzt, um die Auswirkungen der Nanover- stärkung auf die thermomechanischen Eigenschaften des Formgedächtnispolymers zu untersuchen. Weiterhin werden analytische Modelle vorgestellt und zur Berechnung der gesamten effektiven elastischen Konstanten des RVE eingesetzt. Andererseits wird die zelluläre Struktur des SMPF auf der Mesoskala II untersucht. Zwei 3D-Geometrien werden erzeugt und in der FE-Methode verwendet. Desweiteren wird ein rechnerisches Homogenisierungsschema verwendet, um die 3D-RVEs, die auf der unteren Skalenebene erzeugt wurden, in die Analyse einzubeziehen. Zusät- zlich wird ein numerischer Druckversuch mit endlicher Dehnung auf beide Phasen (gummiartig und glasartig) der SMPFs angewandt. In den nächsten Schritten werden alle geschätzten Parameter der mesoskaligen II-Ebene akkumuliert und direkt in die makroskopischen konstitutiven Modelle der Zusammenfassung xxi gummiartigen und glasartigen Phase angewendet. Die thermomechanische Leistung des mit Nanoton gefüllten Formgedächtnis-Polymerschaums wird bei Temperaturen unterhalb und oberhalb von Tg untersucht. Darüber hinaus werden zwei verschiedene thermomechanische Belastungszyklen aufgebracht, um das Rückstellverhalten der SMPF-Nanokomposite zu untersuchen. Schließlich bietet die Multiskalenmodellierung von SMPF ein tiefes Verständnis des Materialverhaltens, insbesondere wenn es von mikro- und nanoskaligen heterogenen Materialstrukturen stammt. Desweiteren wird die Multiskalen-Material-Modellierung (Multiscale Materials Modelling) (MMM) als eines der genauesten Simulationsver- fahren angesehen, welches in der Lage ist, komplexe Strukturen der heterogenen Nanokomposit-Materialien zu modellieren. Zusammenfassung xxii Chapter 1 INTRODUCTION AND REVIEW OF LITERATURE 1.1 Introduction Shape memory materials SMMs are defined as materials that can recover (memorize) their previous form by applying an external stimulus like the thermal or magnetic field. The shape recovery phenomenon (phase transformation) was discovered by Buehler et al. [22] on nickel-titanium NiTi alloys. Many metallic alloys express shape recovery effects when they are subjected to a stimulus (temperature or magnetic field). Solid state phase transformation is used to characterize the shape memory alloy (SMA) such that the starting phase (permanent phase) and the temporary phase are solid structures, Barbarino et al. [14]. Furthermore, SMAs have a wide range of application, especially in microelectromechanical system (MEMS) like thermal and optical sensors, actuators and grippers as well as in Aerospace etc. [160]. The shape recovery effect is not exclusive to metals. Many other materials demonstrate the same response to a stimulus, for instance, shape memory polymers (SMP), shape memory hybrid (SMH) and shape memory get (SMG). Among all these types, SMM and SMP are considered the most important and interesting ones [77]. SMPs were first introduced by Rainer et al. [137] by producing polyethylene (PE) heat-shrinkable films. The mechanism of the shape recovery phenomenon in SMPs 1 1.1. Introduction 2 is related to their intrinsic molecular structure. Thermally induced SMPs show two separated phases. According to the external temperature, SMPs transform from one phase to the other at a specific temperature called transition temperature Tg. Above this temperature, SMPs show the soft phase that stores the permanent shape of the material (Rubbery phase). Whereas, SMPs exhibit a stiffer behaviour at a temperature below Tg and present the glassy phase. Any change in the material from below Tg can be released by heating the material above Tg and return to the original permanent shape. Tg and the elastic properties of every phase are illustrated in figure (1.1) below [49, 98, 105]. Figure 1.1: Elastic modulus variation with respect to the temperature of amorphous polymers [49] In addition to the shape recovery effect, SMP have more advantages [78]: • They can use multi-external stimuli. In addition to heat, SMP has many differ- ent ways to stimulate and recover its original shape (e.g., magnetic field, light and electric field etc.). • SMPs have flexible programming (the process of specifying the permeant shape of the material). Such a process can be done for polymers with different stimuli. 1.1. Introduction 3 • SMPs can be built from a wide range of polymer types. In other words, net- points can be designed and various types of phases can be switched for. • The thermomechanical properties of SMPs can be easily adjusted to meet the requirements by addition of reinforcements or by the synthesis methods. • Their biodegradability and biocompatibility made them very comfortable to produce device items that interface with the human body and bio tissues. Thus, SMPs are very well applicable in the biomedical field. • SMP foams were introduced as a very lightweight material that can carry a high deformation with relatively lower density and higher thermomechanical properties. Such an option made the SMPFs a first choice in the Aerospace field. Despite their superiority in many fields as compared with SMAs, SMPs have low thermomechanical properties. As a result, the application area for SMPs is limited. Thus, nano-filler reinforcements are introduced as sufficient manner to enhance the material without affecting the shape recovery effects. By looking in literature, a number of reinforcing ways are used to improve the material performance. Cao and Jana [23] employed the exfoliated reactive Nano clay particle to improve the shape recovery force of the polyurethanes. Furthermore, Kim et al. [87] prepared a nanocomposite by in situ polymerizations of ethyl methacrylate with intercalated Na-MMT (sodium Montmorillonite). The resulting SMP nanocomposites show good improvement in the mechanical and rheological properties. On the other hand, research groups used carbon nanotubes to increase the mechanical properties of the SMPs. Such nano-filler carbon Nanotubes improved the mechanical and thermal properties as well as the shape recovery force of SPMs dramatically [119, 146, 193]. Many application fields like aerospace and construction prefer materials with low density and higher thermal isolation. Their requirements have motivated the researchers to develop shape memory polymer foams SMPFs which exhibit an impact relaxation property, low density, high thermal isolation, and relatively higher 1.2. Literature survey 4 compressibility [37, 165, 170, 171]. Last but not least, Quadrini et al. [136] prepared nano-clay filled SMPFs that combined the lower density and higher mechanical properties due to nano-reinforcement. Furthermore, a solid-state foaming process, which is suitable for the production of composite foams, is employed to prepare an epoxy foam that depicts enhanced properties with MMT nano-reinforcements. SMPFs Nanocomposites have a heterogeneous structure at individual levels of observation. For instance, at the macroscopic level, the material structure consists of a polymer matrix with irregular and randomly distributed voids. At a lower level of observation, a two-phase epoxy polymer nanocomposite structure is built from a polymer matrix and exfoliated nano-clay layers. It has been generally documented that the core microstructure originates the majority of the macroscopic thermomechanical phenomena. In other words, the main constituents that make up the microstructure, like size, shape, weight fraction are affecting the macroscopic thermomechanical behaviour of the material. Studying the relation between the different scale levels leads to a deep understanding of the material [117]. The main objective of this work is to prepare a full multiscale model that can be used to investigate the impacts of the different scales of material structures on the overall macroscopic SMPFs performance. Furthermore, the effects of upscaling the problem from nanoscale to macroscale level will be examined. 1.2 Literature survey Epoxy shape memory polymer foams depict different heterogeneous structures at a specific observation level. Due to this high heterogeneity, a multiscale modelling method is introduced as a powerful tool to simulate such materials. In the bulk of literature, many research groups employed two- or many-scale continuum mechanical modelling methods to simulate the thermomechanical performance of polymer nanocomposites. Many studies are introduced to establish constitutive models for shape memory polymers. The presented researches targeted the shape recovery phenomena under 1.2. Literature survey 5 small strains. For instance, Liu et al. [102] prepared a constitutive model to investi- gate the thermomechanics of shape recovery of the epoxy polymer. Furthermore, the model is tested within small strains and under uniaxial tension and compression. In like manner, Tobushi et al. [169] developed a thermomechanical constitutive model to explore the thermomechanical properties of SMPs. The model is implemented on the polyurethane series. Thereafter, Westbrook et al. [185] presented a constitutive model that allows to investigate the shape recovery behavior of semicrystalline poly- mers. The estimated model is implemented to illustrate the effects of temperature rates of the actuation strain in a two-way shape memory effect. Sheng et al. [152] predict the overall elastic properties of polymer/clay nanocompos- ites using continuum-based micromechanical models. Since there is no model that represents the material dynamics flow at the atomistic scale, this approach cannot be considered as a full multiscale modelling approach. However, Borodin et al. [19] used a Molecular Dynamic (MD) method to simulate the interfacial polymer shear modulus and include it in the Material Point Method (MPM). By applying MPM, each material point is assigned several material properties, in addition to different properties to particles inside the same background grid element. The study indicated that turning on the attraction between the polymer matrix and nano-inclusions could increase the time-dependent shear modulus. Similarly, a continuum-based micromechanical model was developed by Odegard et al. [122]. The model incorporates the molecular structure of the matrix and nano-particle as well as the interfacial region. The resulting properties of the atomistic scale model are included in the micromechanics model. Such a simulation manner is also used by Yang et al. [189] to characterize the carbon nanotubes size and weakened bonding effect between the polymer matrix and nano-inclusions. Furthermore, the effective elastic stiffness is investigated by using a nanocomposite unit cell. Paliwal and Cherkaoui [126] explored the elastic properties of polymer nanocomposites using an interfacial model. The results are included in an Eshelby type micromechanical scheme. Comparatively, a three-dimensional FE method in parallel with the MD method is used by Choi et al. [31] to establish a multiscale 1.2. Literature survey 6 model which is able to characterize the mechanical properties and the interfacial strength of the polymer nanocomposites. Furthermore, the effect of the particle size on the mechanical behaviour is also investigated. All the above-mentioned studies are focusing on the thermo-elastic properties of the material within a small strain frame, whereas, polymer nanocomposites often undergo finite deformations. Thus, many viscoelastic parameters could not be investigated without testing polymer nanocomposites with finite tension/compression strains. Therefore, a considerable number of studies that developed constitutive models under finite strains have been presented. The earliest model was developed by Tobushi et al. [170], in which, the linear constitutive equations are expanded to nonlinear equations. Consequently, a 3D thermo-viscoelastic model is introduced by Diani et al. [39]. Such a model is established depending on the standard thermo-viscoelasticity theory. At a later time, Chen and Lagoudas [26] expand the constitutive model of Liu et al. [102] to investigate the shape recovery behavior of SMPs undergoing finite strain. A set of nonlinear equations are derived to predict properties of SMPs. At the same instance, a finite deformation constitutive model is introduced by Qi et al. [135]. The model inspects the thermomechanical behavior of thermally induced SMPs that undergo finite deformation. Moreover, the first-order phase transition concept is employed to estimate the change in the material behavior from the elasticity of the rubbery phase to the viscoelasticity of the glassy phase. On the other hand, Nguyen et al. [116] developed a constitutive model that targeted the thermomechanical behavior of amorphous SMPs. The main effects that are included inside the model are structural and stress relaxations. Furthermore, those effects are incorporated in the form of viscoelasticity within the rubbery region and glass transition and viscoelasticity in the glassy region. Later on, a viscoplastic shape recovery constitutive model is presented by Shojaei and Li [154]. The model can capture the temperature-dependent properties within the finite deformation framework. By the same token, a constitutive model with a multi-element is developed by Sweeney et al. [161] to examine large tensile deformations, stress relaxation and stress recovery of SMPs. 1.2. Literature survey 7 The above mentioned finite deformation constitutive models are the most widely employed in modeling of shape memory polymers. For example, a finite element model is implemented by Reese et al. [143] on SMPs at the macro- and microscale level and the thermomechanical response of complex stent structures is investigated. Besides, Volk et al. [180] have reduced the constitutive model proposed by Chen and Lagoudas [26] to one dimensional equations to establish the axial response of SMPs. There, incompressible and isotropic neo-Hookean behaviors are assumed. Over and above that, shape memory cycling tests are modeled by Diani et al. [40]. The thermomechanical behavior of SMPs under large deformation torsional load is explored using the finite element method. On the other hand, the stress-strain response of the shape memory polyurethane is predicted by Pieczyska et al. [133] depending on the Qi et al. [135] constitutive equations. The model is implemented on temperature and strain rate ranges close to room temperature and for affordable strain rates. The power of a multiscale modelling method is fully exploited by Song et al. [155] to simulate not only the material behaviour at finite deformations but also to study the damage progression of nylon 6/clay nanocomposites. A 3D Representative Volume Element (RVE), as well as the MD method, are employed together to build the hierarchical two-scale model. In the same way but without using the MD method, Dai and Mishnaevsky [34] employed multiple disk-shaped nanoplatelets FE unit cell models to analyse the initiation and growth of micro cracks in the nano-clay reinforced polymer composites. Recently, Li et al. [101] presented a modular-based multiscale modelling approach for viscoelasticity of polymer nanocomposites. Furthermore, four modules have been included. Namely, neat polymer toolbox, interphase toolbox, microstructural toolbox, and homogenization toolbox. Both hyperplastic and viscoelastic performance of the polymer nanocomposites are accurately predicted. On the other hand, polymeric foams have acquired a large amount of interest from different researchers. Many modelling studies have been applied to estimate their elastic properties as well as their nonlinear behaviour under finite strain loads. 1.3. Outline 8 Viot et al. [177] has implemented a high-speed camera technique on a dynamic compression device. The evolution of the stress response of polymeric cellular materials has been estimated. SMPF’s are constitutively modelled by Xu and Li [187]. In addition, the effects of classical viscoelasticity and pseudo-plasticity on thermomechanical properties are included in the model. Viot et al. [178] developed a multiscale model able to describe the polymeric foam behaviour under a finite strain compressive test. Furthermore, the main three stages of the global stress-strain response evolution and the local damage by buckling are simulated. Moreover, a new numerical approach with a discrete element method is initiated. Last but not least, Di Prima et al. [38] has generated 3D finite element meshes of the epoxy shape memory polymer foam using CT-images. A hyperelastic model is employed in parallel with the FE method to predict the compression response of the cellular material. In this dissertation, a full multiscale model able to simulate the thermomechanical performance of a nano-clay filled epoxy shape memory polymer foam under finite strains and at atomistic, micro-, meso- and the macro-scale level is presented. In addition, the nonlinear hyperelastic and viscoelastic behaviour of the material is studied as well as its shape recovery behaviour. To the best of the author’s knowledge, no such work has been reported in literature. 1.3 Outline The main objective of this work is to develop a full multiscale model able to simulate the SMPFs thermomechanical performance at different scale levels. Furthermore, a combination of a multiscale FE method, homogenization schemes, RVE concepts, and molecular dynamics method is employed to investigate the material behaviour at different length scales and for different material phases. According to all given motivations, the work schedule of this dissertation can be illustrated as follow: Within chapter 2, the constitutive models of the two different material phases (rubbery and glassy phase) are presented in details. Furthermore, SMPFs prepara- 1.3. Outline 9 tion and the internal structure are illustrated. Finally, the material’s experimental response under finite compressive strain is presented as well. In chapter 3, full molecular dynamics models are proposed to determine the thermomechanical properties of the material nanocomposites at the atomistic scale level. Furthermore, the interfacial strength between the epoxy polymer matrix and the layers is estimated in all possible layer orientations. Moreover, the effects of nano-reinforcements on the thermal expansion coefficient as well as the transition temperature are tested using the MD method. Chapter 4 contains a combination of the FE method and an RVE concept as well as numerical homogenization schemes which are employed to investigate the effects of nano-reinforcements on the thermomechanical properties of the shape memory polymer. Furthermore, two analytical models are employed to estimate the effective elastic constants of the nano-clay filled shape memory polymer. In this chapter, we restrict ourselves to include only two models, namely, the Mori-Tanaka and a Self-Consistent model. They are frequently used for heterogeneous materials for different inclusion shapes. The effects of the cellular structure on the mechanical properties of the SMPFs are illustrated in chapter 5. Two 3D geometries are generated and used in the FE method. In addition, a computational homogenization scheme is used to include the 3D RVE that generated in chapter 4 in the analysis. Moreover, a numerical finite strain compressive test is applied to both phases (rubbery and glassy) of the SMPFs. All the computed parameters from the previous chapters are accumulated and used to explore the SMPFs thermomechanical performance at the macroscopic level in chapter 6. Furthermore, thermomechanical load cycles are applied at the macroscopic level to investigate the shape recovery response. Finally, all the obtained achievements are summarized in chapter 7 in addition to conclusions and future works. 1.3. Outline 10 Chapter 2 THERMOMECHANICAL BEHAVIOUR OF SHAPE MEMORY POLYMER FOAMS 2.1 Material description and experimental observa- tions Epoxy shape memory polymer foam is prepared by Quadrini et al. [136]. Epoxy shape memory polymer is consisting of (3M scotchkote 206N) as a resin with a density of 1.44 gr cm−3. Furthermore, dicyandiamide is used as curative. The mixture takes the powder form and could easily be shaped in tablets. Nano-clay (laviosa delight 43B) was used as filler to enhance the epoxy polymer. The filler is extracted from the naturally occurring montmorillonite MMT. Furthermore, the filler get purification and modification with quaternary ammonium salt to make it suitable for the organic material. The filler particle size was expected to be about 7-9 nm with a real density of 1.6 gr cm−3 [136]. Solid- state foaming processes are employed to produce SMPFs. Such that, the epoxy powder is shaped in tablets using a cylindrical stainless steel mold. In addition, the tablets were fabricated with different MMT weight fractions, namely, 0, 1, 3, 5 and 10 wt %. Afterward, the tablets are inserted into a muffle and get foamed at 320 ◦C for 8 min.. Serial experimental tests are performed on SMPFs samples, namely, compression, flexure, stress relaxation, and indentation. 11 2.1. Material description and experimental observations 12 Compression test are performed on a cylindrical foam sample with 10 mm height and 20 mm diameter. Furthermore, 5 mm/min deformation rate is used with a maximum strain of about 80 % at room temperature. The results of different MMT wt % are presented in figure (2.1) below [136]. Figure 2.1: Compression curves for all the foams for 30 ◦C [136] In the same manner, the compression test is performed at a high temperature (130 ◦C) to investigate the material performance with the rubbery phase. Figure (2.2) shows the stress-strain curves for various MMT weight fractions [157]. Moreover, a shape recovery test is also performed on the foam sample with MMT weight fractions of 5 % the material showed a 97 % height recovery as depicted in figure (2.3) and figure (2.4). The final cellar structure of the material is shown in figure (2.5). Due to the physical solid-state foaming at high temperature, the microscopic structure of the material is expected to be fully exfoliated silicate layers distributed randomly through the epoxy matrix. 2.1. Material description and experimental observations 13 Figure 2.2: Compression curves for all the foams for 130 ◦C [157] Figure 2.3: Stress as a function of time and strain during the packing and cooling phases of a thermo-mechanical cycle (MMT content: 1 wt%) [136] 2.2. Multiscale modeling approaches 14 Figure 2.4: Stress Relaxation test [136] Figure 2.5: Left: Surface morphology of the foam structure adopted from [157], right: the expected exfoliated nanostructure adopted from [115] 2.2 Multiscale modeling approaches It is clear from the previous section that the SMPFs have different heterogeneous structures at different scale levels. At a specific level, the material has a cellular structure, in which, many voids with different sizes are distributed randomly through the epoxy polymer matrix. Furthermore, open cell foams have a rough internal cell surface and a very high irregularity. On the other hand and at a lower level, the material structure is built from two phases, namely, epoxy polymer matrix and fully exfoliated silica layers. The filler’s layers have different orientations as shown in figure (2.5). However, preparing one numerical model which is able to simulate the material thermomechanical behavior and include all the effects of the lower scale material structure which is a very hard task. Alternatively, a multiscale modeling approach is 2.2. Multiscale modeling approaches 15 Figure 2.6: Heterogenous structure of shape memory polymer foams at different scales employed to simulate SMPF at various scale levels. A multiscale modeling method is performed by dividing the problem into four separated moods according to the scale level as shown in figure (2.7). Furthermore, numerical and computation homogenization schemes are used to link and upscale the various scale level model. Moreover, in each level, a numerical model is applied and the effects of the material structure at that level on the material behaviour are studied in detail and are used in the next length scale level. Starting from the atomistic scale level, molecular dynamics method MD is used to estimate the elastic parameters of the material as well as the interaction strength and traction-separation curves between the epoxy matrix and the nano-clay fillets. The obtained parameters and results are used in the next upper scale level mesoscale I. The effects of the MMT silicate layers on the overall effective elastic parameters are investigated by using the finite element method and RVE concepts. Furthermore, the RVEs are tested with finite deformation compressive loads. At the mesoscale II level, a computational homogenization scheme is used to include the effects of the nano-reinforcements on the mechanical performance of the SMPFs. Furthermore, the effects of a cellular structure and the density reduction on the physical properties of the material are explored. Finally, all the estimated parameters are accumulated and applied to the constitutive models of the material at the macroscopic level. Such that, finite strain numerical compressive tests are performed on both rubbery and glassy phases. Last but not 2.3. Constitutive models 16 Atomistic scale 100 Å Mesoscale I 50 nm Mesoscale II 10 mm Macroscale 10 cm Length sclae T im e sc la e -Elastic constants -Adhesiv e stre ngth (Wad h ) -Glass transiti on (Tg ) -Expansion coefficien t (α) -Physica l properit ies (E ∗ , G ∗ K ∗ , v ) -σ − ε curves -Lamé constants (λ, µ) -Initia l hardening muldulus (Cr) -Soften ing slope (h) -Stea dy state value (S0 ) -Prefe rred state value (Sss ) Figure 2.7: Schematic of a multiscale modeling approach for SMPFs nanocomposites least, the shape recovery behavior of the material is tested using thermomechanical load cycles. In the next chapters, full detailed descriptions of every scale level are presented. 2.3 Constitutive models The thermomechanical response of the shape memory polymers is highly temperature dependent. At room temperature, SMPs show a stiff response. Such SMP can become much softer and behave like ductile rubber at higher temperatures. The transition temperature Tg separates these two different behaviours, figure (2.8), such that, above Tg, the material behaves like a rubbery phase with a nonlinear hyperelastic response. Whereas, below Tg glassy phase become the main phase in SMPs [102]. According to the above-mentioned observations, three constitutive models are needed to simulate the material. First is a constitutive model that is able to capture the nonlinear 2.3. Constitutive models 17 hyperelastic behavior of the rubbery phase at a temperature higher than Tg. Second is a model that represents the nonlinear viscoelastic response of the glassy phase which appears at a temperature below Tg. Finally, a constitutive model controls the transition between the two phases and captures the shape recovery effects. Figure 2.8: Storage modulus, loss modulus and tan delta of the shape memory poly- mer [102] In this approach, the constitutive model developed by Arruda and Boyce [11] is used to describe the nonlinear hyperelastic behaviour of the rubbery phase. On the other hand, the constitutive model developed by Boyce et al. [20] is employed to define the nonlinear viscoelastic response of the glassy phase. Qi et al. [135] have prepared a constitutive model that is able to simulate the phase transition and shape recovery effects of the SMP under finite deformation. Such a model will be used in the following chapters. 2.3.1 Rubbery phase An eight-chain model [11] can capture the rubber-like phase hyperelastic response, starting from an equation relating the Cauchy stress tensor T r with the deformation 2.3. Constitutive models 18 gradient tensor F T r = Cr √ N λchain L −1λchain√ N [ B − λchain 2I ] (2.1) where Cr is the initial hardening modulus, it is a material property defining the strain hardening. N represents the number of rigid molecular units between entan- glements. B = FF T (2.2) λchain = [ 1 3 trB ]1/2 (2.3) λchain is the stretch on every chain in the network. Furthermore, the Langevin function L is given by: L (β) = coth (β)− 1 β (2.4) 2.3.2 Glassy phase Capturing the nonlinear viscoelastic behavior of the glassy phase required three con- stitutive elements as shown in figure (2.9). The linear spring element simulates the initial elastic response of the material. Furthermore, the viscoelastic dashpot used to capture the rate dependent yield. This element is connected in parallel with a nonlinear rubber elasticity spring element which develops the back stress. Starting with the multiplicative decomposition of the deformation gradient[20] F = F eF p (2.5) F p represents the deformation gradient of the relaxed configuration which can be estimated by elastically unloading to a stress-free state by F e−1 . In addition the elastic deformation gradient can be also decomposed into two tensors, stretch, and 2.3. Constitutive models 19 Figure 2.9: Spring-dashpot elements representation of glassy phase constitutive model rotation F e = V eRe (2.6) Moreover, the velocity gradient L can be represented as L = Ḟ F−1 = D +W = Le + F eLpF e−1 (2.7) where D is the rate of stretching and W is the spin. However, the velocity gradient of the relaxed configuration Lp is illustrated as: [Boyce 2001] Lp = Ḟ p F−1 = Dp +W p (2.8) where Dp is the shape change rate of the relaxed configuration, W p is the spin and it will be taken as W p = 0. The Cauchy stress tensor T g can be related with strain as T g = 1 J ϑe [lnV e] (2.9) 2.3. Constitutive models 20 where V e is the hencky strain and J = detV e ϑe = 2µgI + λgI ⊗ I (2.10) where µg and λg are lame constants, I is the fourth order identity tensor and I is the second order identity tensor. Now the rate of shape change of the relaxed configuration can be constitutively prescribed as Dp = γ̇pN (2.11) where γ̇p is the plastic shear strain rate that defined as γ̇p = γ̇0 exp [ −∆G kT ( 1− τ s )] (2.12) where γ̇0 is the pre-exponential factor, ∆G is the zero stress level activation energy, k is the Boltzmann’s factor and T is the absolute temperature. N is introduced as a normalized tensor aligned with the deviatoric driving stress state N = 1√ 2τ T ∗′ (2.13) The deviatoric portion of the driving stress state T ∗′ is presented as T ∗′ = ReT [ T − 1 J F eTNF eT ] Re (2.14) where, TN is defined as the network stress that captures the effects of orientation- induced strain hardening [20]. The effective equivalent shear stress τ can be estimated in the following form: τ = [ 1 2 T ∗′ • T ∗′ ]1/2 (2.15) In addition, the resistance s is employed to evolve with the plastic strain from an 2.3. Constitutive models 21 initial value s0 to the steady-state value by ṡ = h [ 1− s sss ] γ̇p (2.16) where, h is the softening slope. The last equation describes the strain softening, while s is the athermal shear stress and sss is the preferred state. 2.3.3 Shape recovery effect Such a constitutive model is developed by Qi et al. [135]. It supposes that the total strain is divided into the two phases that made the whole material structure (rubbery and glassy) according to their volume fraction fr and fg respectively. The volume fraction of the two phases should satisfy fr + fg = 1 (2.17) such that the total Cauchy stress tensor is determined as T = frT r + fgT g (2.18) The volume fraction of each phase is changed from 0 to 1 according to the tem- perature variation. The volume fraction as a function of temperature is defined as fr = 1 1 + exp [−(T − Tr)/A] (2.19) fg = 1− 1 1 + exp [−(T − Tr)/A] (2.20) where, A is the width of the phase transition zone, Tr is a reference temperature and it is close to Tg. It is noteworthy that during cooling, the glassy phase will be grown. The recently formed glassy phase does not take over the deformation of the rubbery phase but it behaves as an undeformed configuration [135]. In addition, the 2.3. Constitutive models 22 increment deformation gradient of the glassy phase is defined as ∆F g = F n+1F n−1 , if ∆T ̸= 0 I, if ∆T = 0 (2.21) where, F n and F n+1 are the deformations at the increment n and n+1 respectively. However, the total deformation gradient acting on the glassy phase is F n+1 g = ∆F n+1 g F n g (2.22) Last but not least, all the mentioned constitutive models are implemented in ABAQUS/explicit user subroutine VUMAT and applied on the micro- and macro- scopic scale level in the next chapters. Chapter 3 MOLECULAR DYNAMICS SIMULATION OF NANO-CLAY FILLED SHAPE MEMORY POLYMERS 3.1 Introduction After the revolution of the nanotechnology, many nanomaterials are used as inclu- sions to reinforce many types of polymers. Such enhancements increased the interest in polymer nano-composites due to their improved mechanical properties. Recently, the range of application of polymer nano-composites has been extended in many fields, especially in polymer industry and in biomedical applications. The addition of micro-scale particles to shape memory polymers reduces the shape recovery behavior of the material. Thus, developing of shape memory polymer nano-composites is con- sidered as a powerful solution to improve the thermomechanical properties without affecting the shape recovery nature. On the other hand, new problems have arisen. For instance, how to characterize the deboning process between the polymer matrix and the nano-inclusions? Such a process highly affects the mechanical properties of nano-composite materials. More- over, what are the effects of the orientations of the nano-inclusions on the overall behavior of the material in a tension or compression test? Due to the tiny size of the nano-inclusions, an inspection of the dynamics of the material experimentally is really an arduous task. Therefore, molecular dynamics has come to lights as an effective 23 3.1. Introduction 24 simulation method to study atomistic scale dynamic material phenomena which are hard to inspect experimentally. Molecular dynamics has the abilities to simulate material behavior at the atomistic scale. It is a valuable tool to estimate many thermomechanical properties of a mate- rial, which are difficult to measure experimentally. In this chapter, we are focusing on the elastic constants, transition temperature, and the coefficient of thermal expansion. Furthermore, the interfacial strength at the interface zone between the matrix and the nano-inclusions can also be determined using molecular dynamics. The resulting parameters are used as base parameters for the higher scale continuum thermome- chanical model. Recently, many attempts have been reported in literature on using molecular dynam- ics to simulate epoxy/clay nano-composites at the atomistic scale. Zhang et al. [195] employed molecular dynamics together with the inverse gas chromatography (ICG) to analyze the effects of the polarity of the organic modifier and polymer matrix on morphology and interfacial strength of the polymer/clay nano-composites. Nguyen et al. [118] presented a molecular approach to investigate the morphology and ther- momechanical properties of epoxy/clay nano-composites. The study has been applied to intercalated and exfoliated structures of epoxy/clay nano-composites. Chen et al. [27] introduced a full study on the debonding behavior between the matrix and the nano-clay. The traction separation curves of both the clay and matrix interface are obtained. Moreover, Xu et al. [188] predicted the overall compressive moduli on Ny- lon 6/MMT nano-composites using the MD method. The study contained fully and partially exfoliated clays structures. Likewise, the effects of clay cluster alignments on the effective compressive moduli were investigated. In a similar manner, Tanaka and Goettler [162] and Toth et al. [174] concentrated on the prediction of the binding energy among Nylon 6, organic modifiers and clay platelet. They used molecular dynamics simulations to conclude that there is a relationship between the resulting binding energies and the fracture energy between the Nylon 6 matrix and the clay layers. In the same fashion, a further study was done by Kim and Kim [86]. In the mentioned 3.2. Material description 25 work, the moisture diffusion characteristics and mechanical properties of epoxy/clay nano-composites were investigated using molecular dynamics simulation as well as experiments. Additionally, the tendency of the mechanical performance improve- ments was successfully predicted. An additional study, which focuses on the effects of nano-inclusions on the glass transition temperature and the coefficient of thermal expansion, was presented by Choi et al. [30]. They applied molecular dynamics on epoxy nano-composites containing spherical nano-sized SiC. Clear improvements in the studied parameters have been observed. After a deep literature survey, the ma- jority of the presented studies are focusing on the interfacial strength and the fracture energy between the polymer matrix and the nano-clay in the normal direction (the debonding is perpendicular to the silicate layer surface). To our best knowledge, the presented study is the first work that studies the interfacial strength at the atomistic scale not only in normal direction but also in the shear direction and all the possible orientations of the silicate layer alignment. Furthermore, full molecular models ap- proaches are presented to estimate the thermomechanical properties of the material nano-composites. 3.2 Material description The first step of a molecular dynamics simulation is to find out the molecular structure of all the chemical compounds that are present inside the materials system. As a nano-composite material, nano-filled shape memory polymer foams are built from three main phases, namely, epoxy polymer (matrix), silicate layers (inclusions) and organic modifier (interphase). Within this section, the molecular structures of the three mentioned phases will be discussed in detail. The epoxy polymer (matrix) phase is generated from the polymerization process, in which, DGEB A resin is cross-linked by addition of Dicyandiamide as a harder and curing (cross-linking) agent. The molecular structure of DGEB A consists of a group of molecular chains. Every chain involves a number of repeated monomers with two epoxy rings at the ends as shown in figure (3.1) [35, 54, 124]. Furthermore, the molecular structure of the Dicyandiamide, which is used as a 3.2. Material description 26 Figure 3.1: DGEB A chains molecular structure curing agent in the material, is shown in figure (3.2). It is commonly used in the heat- curing epoxy resin for coating and adhesive. In addition, it cures epoxies speedily at elevated temperatures (120-180 ◦C) and has a high purity as well as low reactivity at a temperature below 100 ◦C. The conventional weight fracture of Dicyandiamide that is used to formulate epoxy resin is 3-10 [67]. Figure 3.2: The molecular structure of the Dicyandiamide The curing mechanism of DGEB A resin by amine is described by Gilbert [60]. The equations of the curing chemical reactions are illustrated in figure (3.3). Such a reaction produces a three dimensional cross-linked epoxy network. Generally, a highly cross-linked network leads to a low deformation capability and a relatively high heat resistance [88]. Figure 3.3: The reactions that occur between resin and curative agent As mentioned in chapter 2, nano-clay (Laviosa Dellite 43B) was used to reinforce 3.2. Material description 27 the epoxy polymer system. It is derived from montmorillonite (MMT) which occurs naturally. Due to its highly charged layers, MMT needs some purification and modi- fication processes to reduce its high charge and to make it more compatible with the epoxy polymer. The molecular structure of the Na+MMT consists of 1 nm thick lay- ers. Every layer is built from two tetrahedral sheets of silicate separated by one sheet of an edge-shaped octahedral with alumina or magnesia, figure (3.4). The layers are stacking and that causes a gallery or an interlayer which is nothing but an orderly Van der Waals gap between the silicate layers. The negative charges, that are generated from the isomorphic substitution inside the layer, are commonly counterbalanced by hydrated Na+ cations that are existing in the interlayer [58]. Figure 3.4: Molecular structure of the silicate layers, adopted from [58] The existence of Alkyl Na+ cations between the layers makes MMT immiscible with polymers and organic materials. Thus, exchanging this cation with an organic ion like ammonium is very desired for two reasons: first, the organic ion with its high molecular weight increases the gap between the MMT layers through the intercalation and relatively increase the internal surface area. Second, the inserted alkyl ammonium cation will increase the compatibility of MMT with polymers because it has the ability to react with the polymer or entangle with the polymer chains. Absolutely, that will improve the adhesion strength between the epoxy matrix and the silicate layers [8, 130, 141]. Moreover, separating the silicate layers by intercalation give the polymer chains more chances to exfoliate the layers from each other during the polymerization process [28, 95]. The expected final structure of the material shows the polymer matrix with randomly distributed and oriented silicate layers as shown in figure (3.5) 3.3. Molecular Dynamics (MD) Method 28 below. Figure 3.5: Structure of the polymer/nano-clay composites with fully exfoliated lay- ered silicate, adopted from [29, 115] Figure (3.6) expresses the intercalation and exfoliation during the modification and polymerization processes. The organic modifier (om) that is used to modify the MMT is quaternary ammonium salt (dimethyl benzyl hydrogenated tallow ammonium). The molecular structure of the organic modifier is illustrated in figure (3.7). The main benefits of using such om are their relatively high molecular weight and the existence of benzene rings which can play a big role in the intercalation and exfoliation processes. Furthermore, om establishes an entanglement with the polymer chains and improves the interface strength between the polymer matrix and the silicate layers. 3.3 Molecular Dynamics (MD) Method Consider a 2N-dimensional phase with N particles with masses {m1,m2, .....mN}. The system is represented by the coordinates of particles {x1,x2, .....xN} along with the momenta p := {p1,p2, .....pN} T . The dynamic of such system is governed by the Hamiltonian equation as [94], H (x,p) = N∑ i=1 p2 i 2mi + U (xi) (3.1) where U (x) is the potential that characterizes the interaction between the particles. 3.3. Molecular Dynamics (MD) Method 29 Figure 3.6: Representation of the intercalation and exfoliation of the silicate layers during the modification of the nano-clay and polymerization process Figure 3.7: The molecular structure of the organic modifier (om) that is used in the purification and modification of the nano-clay Furthermore, the Hamiltonian equation of motion can be described as, dxi dt = ∂H ∂pi dpi dt = ∂H ∂xi (3.2) The mentioned equations of motion (3.1) and (3.2) are equivalent to Newton’s equa- tions of motion that are illustrated by Newton’s second law [65], miẍi = F i (x, t) (3.3) 3.3. Molecular Dynamics (MD) Method 30 where F i is the force acting on the particle i and is defined as, F i = −▽xi U (3.4) 3.3.1 Statistical ensembles 3.3.1.1 NVE ensemble The NVE ensemble is known also as the microcanonical ensemble. It is obtained by solving Newton’s equations of motion with keeping temperature and pressure varying along the trajectory. In this ensemble, the number of particles, volume, and total energy stay constant along the time trajectory. Introducing the total energy of the system which results from the summation of the kinetic energy Ekin and potential energy U is given as [65, 76, 129, 139, 140] Etot = Ekin + U (3.5) Here, the kinetic energy can be defined as [65], Ekin = 1 2 N∑ i=1 pT i pi mi (3.6) Let α be a dynamic variable, then the time average can be estimated due to the ergodic hypothesis [10] as, < α >t= lim t→∞ 1 t ∫ t 0 α (τ) dτ (3.7) Here, α can be replaced by any function of coordinates and momenta of the parti- cles system. Using the equipartition theorem, the thermodynamic temperature T is defined by averaging the kinetic energy at a specific time t as [65], < Ekin >t= NfkBT 2 (3.8) where, kB is the Boltzmann’s constant and Nf is the physical system’s degree of 3.3. Molecular Dynamics (MD) Method 31 freedom. The temperature can then be illustrated as [65], T = 2 NfkB Ekin (3.9) In particular, the thermodynamic pressure P could be calculated. Firstly, let’s introduce a relation between the thermodynamic pressure and the internal virial ⟨W ⟩t by using the classical virial theorem [5], such that [65], PV = 1 3 NfkBT + ⟨W ⟩t (3.10) The variation of the potential energy with respect to volume is also related to the internal virial W by W V = − dU dV . Then we can defined the pressure in the form [65] P = 2 3V Ekin + 1 V W (3.11) Using the above equations, we can now implement Molecular Dynamics simula- tions, in which, N, V, and E are kept constant. Furthermore, in the microcanonical ensemble, the preferred temperature cannot be reached because there is no energy flow aided by the temperature control methods. Thus, the above ensemble is not recommended for equilibration and minimization processes. 3.3.1.2 NPT ensemble It is also known as the isothermal-isobaric ensemble. In this ensemble, the number of particles, pressure, and temperature are constants. In addition, it provides a full control over temperature and pressure. Furthermore, the volume of the unit cell has the ability to change during the simulation. Therefore, the pressure is adjusted by adjusting the unit cell vectors. In order to establish relations that control the thermodynamic quantities within the NPT ensemble, let’s introduce the Lagrangian in the form [65, 120, 121, 127] L (x,v) = 1 2 N∑ i=1 miv T i vi − U (xi,x2, ....xN) (3.12) 3.3. Molecular Dynamics (MD) Method 32 where, xi represents the coordinates of the particle’s position and vi are the velocities of the particles. In addition, the equation of motion is given by the Euler-Lagrangian equation [65], ▽xL (x,v) = d dt ▽vL (x,v) (3.13) The momenta in this space are pxi = ▽vi L (x,v) = mivi (3.14) by inserting equation (3.14) into the Hamiltonian, we obtain [65] H (x,px) = N∑ i=1 pT xi pxi mi − L ( x1, ...xN , pxi m1 , ....., , pxN mN ) = 1 2 N∑ i=1 pT xi pxi mi + U (x1, ...xN) (3.15) The main target is to formulate an isothermal-isobaric (NPT) ensemble with the possibility of controlling pressure P and temperature T. However, coupling the iso- lated physical system to a faked external pressure and heat is needed to build such ensemble. Therefore, additional degrees of freedom for the total coordinate system in space are developed to make the volume and the shape of the system’s domain and thus the pressure controllable [127, 128]. Similarly, we need more dynamical variables to control the temperature of the system. Let’s assume a unit cell in space R3 at the initial undeformed configuration with standard orthonormal basis e1, e2 and e3 To find dynamical variables that control the shape and volume of the unit cell, we use a 3× 3 time-dependent cell transformation matrix h̃ = [a1, a2, a3] such that [65], xi = h̃s̃i (3.16) where s̃i is the scaled coordinate vector of the particle i. Now, the dynamical variables are necessary to control the temperature of the system. Thus, the time t is rescaled 3.3. Molecular Dynamics (MD) Method 33 to t̃ by [65] t ( t̃ ) = ∫ t 0 dτ γ (τ) (3.17) which leads to [65] dt̃ = γ ( t̃ ) dt (3.18) In the same manner, the velocities of the particles can be scaled as [65] vi (t) = γ ( t̃ ) h̃ ( t̃ ) ˙̃si ( t̃ ) (3.19) After scaling of the variables t and v, they are divided into two groups, namely, virtual variables labeled by tilde and real variables given in terms of real time designated without tilde. In addition, further fictitious potentials of the dynamic variables can be defined, thus, the equations (3.9) and (3.10) can be rewritten as [65] UP ( h̃ ) = Pdet ( h̃ ) (3.20) and UT (lnγ) = NfkBT lnγ (3.21) where, P is the external pressure and T represents the targeted temperature. Subsequently, the Parrinello-Rahman-Nose Lagrangian of the NPT ensemble for one degree of freedom γ is defined as [65] LNPT ( s̃, h̃, γ, ˙̃s, ˙̃h, γ̃ ) = 1 2 N∑ i=1 miγ 2 ˙̃sTi ˙̃hT ˙̃h ˙̃si + 1 2 MPγ 2tr ( ˙̃hT ˙̃h ) + 1 2 MT γ̇ 2 − U ( h̃s̃, h̃ ) − Pdet ( h̃ ) −NfkBT lnγ (3.22) The MP mentioned above is the inertia parameter. It is used to control the motion’s time-scale of cell h̃. MT is a similar parameter with respect to the temperature. 3.3. Molecular Dynamics (MD) Method 34 Corresponding to equation (3.14), the conjugated momenta could be introduced as [65]: ps̃i = ▽ ˙̃si LNPT = miγ 2G̃ ˙̃si (3.23) ph̃ = ▽ ˙̃ h LNPT = γ2MP h̃ (3.24) pγ = ▽γ̇LNPT =MT γ̇ (3.25) where G̃ is the metric tensor G̃ = h̃T h̃. Furthermore, submitting the above momenta in equation (3.15), the Parrinello-Rahman-Nose Hamiltonian gives [65]: H̃NPT ( s̃, h̃, γ,ps̃, ph̃, pγ ) = 1 2 N∑ i=1 pT s̃i G̃−1ps̃i miγ2 + 1 2 tr ( pT h̃ ph̃ ) γ2MP + 1 2 p2γ MT (3.26) Similarly, equation (3.2) is used together with equation (3.26) to formulate the equa- tions of motion as follows [65], ds̃i dt̃ = ▽ps̃i H̃NPT = G̃−1ps̃i miγ2 dps̃i dt̃ = −▽s̃iH̃NPT = −▽s̃iŨ ( s̃, h̃ ) dh̃ dt̃ = ▽ph̃ H̃NPT = ph̃ γ2MP dph̃ dt̃ = −▽h̃H̃NPT = −▽h̃ ( 1 2 N∑ i=1 pT s̃ G̃ −1ps̃ miγ2 ) − ▽h̃Ũ ( s̃, h̃ ) − ▽h̃UP ( h̃ ) dγ dt̃ = ▽pγH̃NPT = pγ MT dpγ dt̃ = −▽γH̃NPT = 1 γ  N∑ i=1 pT s̃i G̃−1ps̃i miγ2 + tr ( pT h̃ ph̃ ) miγ2 +NfkBT  (3.27) 3.3. Molecular Dynamics (MD) Method 35 where Ũ is represented as [65] Ũ ( s̃, h̃ ) = U ( h̃s̃, h̃ ) (3.28) Up to the present time, we have got the Hamiltonian (3.26) coupled with the equa- tions of motion (3.27). Although the resulting equations could be solved numerically, still the equations are represented in term of the virtual time t̃. That will lead to a variable time step in the time discretization scheme. To overcome such difficulties, the equations of motion (3.27) should be transformed into the real-time. Thus, new variables relative to real-time are introduced as [65] si (t) = s̃i ( t̃ ) (3.29) h (t) = h̃ ( t̃ ) (3.30) taking the derivative with time give [65], dsi dt (t) = γ ( ˜̃t ) ds̃i dt̃ ( t̃ ) (3.31) dh dt (t) = γ ( t̃ ) dh̃ dt̃ ( t̃ ) (3.32) Consequently, the associated momenta can also transformed as [65], psi = G̃−1ps̃i γ (3.33) ph = ph̃ γ (3.34) 3.3. Molecular Dynamics (MD) Method 36 and their real time derivative as [65] dpsi dt = d ( G̃−1ps̃i 1 γ ) dt̃ = G̃−1 ( t̃ )(dG̃ dt̃ ( t̃ ) G̃−1 ( t̃ ) + dps̃i dt̃ ( t̃ ) − ps̃i ( t̃ ) γ ( t̃ ) ) (3.35) dph dt = γ ( t̃ ) d(ph̃ 1 γ ) dt̃ ( t̃ ) = dph̃ dt̃ ( t̃ ) − ph̃ ( t̃ ) γ ( t̃ ) dγ dt̃ ( t̃ ) (3.36) To express the equations of motion in simple terms, the logarithm of γ is taken as [65] η (t) = lnγ ( t̃ ) (3.37) pη (t) = pγ ( t̃ ) (3.38) thus dη dt (t) = γ ( t̃ ) d dt̃ lnγ ( t̃ ) = dγ dt̃ ( t̃ ) (3.39) dpη dt (t) = dpγ dt̃ ( t̃ ) (3.40) Starting from the equations of motion (3.27) along with the equations (3.29, 3.33, 3.3. Molecular Dynamics (MD) Method 37 and 3.37) the equations of motion in term of real time are expressed as [65] ṡi = psi mi ḣ = ph MT η̇ = pη MT ṗsi = −h−1▽xi U −G−1Ġpsi − pη MT psi ṗh = −▽xi UsTi − ▽hU + N∑ i=1 mihṡiṡ T i − h−TPdet (h)− pη MT ph ṗη = N∑ i=1 pT si Gpsi mi + tr ( pThph ) MP −NfkBT (3.41) Moreover, the transformed Hamiltonian could also be obtained in term of real time in the same manner [65] HNPT (s, h, η,ps, ph, pη) = 1 2 N∑ i=1 pT si Gpsi mi + 1 2 tr ( pThph ) MP + 1 2 p2η MT + U (hs, h) + Pdet (h) +NfkBTη (3.42) The hamiltainian equatin of NPT remains constant over time. Furthermore, the kinetic energy [65] Ekin = 1 2 N∑ i=1 pT si Gpsi mi (3.43) and the potential energy [65] Epot = U (hs, h) (3.44) are not conserved over time. Similar to equation (3.6) the thermodynamic tempera- ture is represented as [65] T = 2Ekin NfkB (3.45) 3.4. Model description and calculations 38 In addition, the thermodynamic pressure is related to the fluctuations of the potential energy U with respect to volume V as [65] P = 2 3V Ekin + 1 V W = 2 3V Ekin + dU dV (3.46) using the chain rule in parallel with the relation V = deth, P is defined as [65] P = 1 3 tr (Πint) (3.47) where, Πint is the internal stress tensor, such that [65] Pint = 1 det (h) N∑ i=1 mihsis T i h T +Πpot int (3.48) Πpot int = 1 det (h) Fhh T (3.49) Fh = N∑ i=1 F is T i (3.50) 3.4 Model description and calculations 3.4.1 Coarse-grain parametrization and Force-Field Introducing a huge number of atoms in the simulation cell is important to get more ac- curate results but this will be computationally expensive. For this reason, it is helpful to reduce the number of degrees of freedom by converting two or three bonded atoms to one superatom. Of course, that should be done in parallel with the availability of United atoms force-field. Figure 3.8 shows the mapping scheme of the United atoms on our model monomers. The OPLS (Optimized Potentials for Liquid Simulations) force fields, developed by [84], are applied in the presented work. They are comfortable with organic ma- terials like polymers. Furthermore, the availability of OPLS United Atoms makes them computationally not expensive. Many physical quantities that can be approx- 3.4. Model description and calculations 39 Figure 3.8: Mapping scheme of the united atoms imated using the Molecular Dynamics method are controlled by the total energy of the atoms inside the simulation box. Such energy is a summation of force-fields that are distributed among the atoms inside the simulation cell. Generally, most of the force-fields have four main forms namely: bond stretching, angle bending, torsional energetics and non-bonded interactions [42, 66, 84, 108, 184]. The representations of the mentioned potentials are illustrated below. The potential of the bond stretching between two bonded atoms is given by equation (3.51) Ebond = ∑ bonds Kr(r − r0) 2 (3.51) where r0 is the zero-force bond distance and Kr is a constant with units of (energy/distance2). The energetics for angle bending among three bonded atoms is Eangle = ∑ angles Kθ(θ − θ0) 2 (3.52) where θ0 is the zero-force value of the angle and Kθ is a constant with units of (energy/radian2). The torsional energy can be written as Etorsion = ∑ torsions K1 2 [1 + cos(ϕ)] + K2 2 [1− cos(2ϕ)] + K3 2 [1 + cos(3ϕ)] + K4 2 [1− cos(4ϕ)] (3.53) 3.4. Model description and calculations 40 where K1, K2, K3 and K4 are constants with energy units. The last intra- molecular potential is the non-bonded interaction between all pairs of atoms (i, j). It is represented by the Coulomb plus Lennard-Jones terms in equation (3.54) Enonbond = qiqj ϵr + 4ϵ [(σ r )12 − (σ r )6] (3.54) where qi and qj are the charge on the atoms i, j respectively, ϵ is the depth of the potential at the minimum r, and σ is defined as the zero-crossing distance for the potential. On the other hand, ClayFF force-fields which were developed by Cygan et al. [33] are employed to define the interactions between atoms inside the silicate layers. Similar to OPLSA, the total energy consists of Coulombic Interactions, Van der Waals interactions, and the bonded interactions Etot = Ecoul + EV DW + Ebondedstretch + Eanglebend (3.55) The energy of interaction in coulombic term is illustrated in equation (3.56). It is inversely proportional to the distance between two atoms rij: Ecoul = e2 4πϵ ∑ i ̸=j qiqj rij (3.56) where, qiqj are the partial charge and they are estimated from the quantum me- chanics equations, e is the charge of the electron and ϵ0 represents the vacuum dielec- tric permittivity. The second term of the nonbonded energy (frequently referred to as the Van der Waals energy) is defined by Lennard-Jones (12-6) function. It involves the short- range repulsion which results from the increase in attractive energy when two atoms are approaching each other. The attractive dispersion energy represented as, EV DW = ∑ i ̸=j D0,ij [( R0,ij rij )12 − 2 ( R0,ij rij )6 ] (3.57) where D0,ij and R0,ij are parameters and they can be determined by the fitting 3.4. Model description and calculations 41 of the model to the physical properties. Furthermore, the geometric mean rule is applied to estimate the energy parameter D0 such that R0,ij = 1 2 (R0,i +R0,j) (3.58) D0,ij = √ D0,iD0,j (3.59) As reported by literature and by experimental observations, there exist no bonds between the silicate layers and the polymer matrix. There are only nonbonded inter- actions between them. That was the main reason for selecting ClayFF in our model for silicate layer atoms. According to [33], the mean of ionic (non bonded) defini- tion of the metal-oxygen interactions is used as a base for the ClayFF. Moreover, all atoms serve as point charge and have complete translation freedom in this force field framework. ClayFF focuses on the nonbonded interactions between the atoms. Such an approach enriches the molecular dynamics model and leads to more acceptable re- sults. In addition, the availability of ClayFF in LAMMPS makes it very comfortable and smoothly applicable in many molecular dynamics models. 3.4.2 Model preparation The main and most important step in MD simulations is to prepare the simulation cell which can be defined as a box that contains all the atoms inside. Furthermore, every atom must align in a specific position with a minimum energy relative to the surrounding atoms. Preparing such a simulation box is not so easy. Unlike metals, the atomic structures of polymers consist of a huge number of molecules with different lengths. These molecules are entangled with each other and are distributed randomly through space. As a Molecular Dynamics solver, LAMMPS is unable to create such complex molecular structures. As a result, other preprocessing tools are needed to build such a model. Moltemplate tool [173] is introduced as a general molecule and force field database system builder for LAMMPS. Furthermore, it provides a powerful possibility to build 3.4. Model description and calculations 42 a force-field database system from the popular force-fields like AMBER, GAFF, and OPLSA-AA. Moltemplate has some limitations. It cannot perform any calculations for the atom positions. It needs also another tool for visualization like vmd [168]. Within this work, the DEBA resins chains are created and distributed inside the simulation box using Moltemplate tool in parallel with the self-avoiding random walk algorithm which was developed by Theodorou and Suter [166]. To create every chain, a non-occupied position Rinit was selected randomly inside the simulation box. Then the first monomer of the chain is created at the base position R0 = Rinit. In a next step, the created monomer is rotated around x, y and z-axes with randomly selected angle θi to calculate the new monomer tip position Rn, such that, −90◦ ≤ θi ≤ 90◦ i=x, y, z (3.60) After checking the occupancy of the new position Rn, it is considered as a base position for the next monomer. The above steps are repeated for every monomer until reaching the desired molecular weight of the chain. The full algorithm is presented by figure (3.9). As mentioned before, Moltemplate cannot perform any calculations. Therefore, the algorithm was applied in a FORTRAN code separately. After creating the resin chains and the curative molecular together inside a (300× 300× 300)Å simulation cell, the cross-linking process is required to bond the chains together and build the final epoxy network. The carbon atoms at every resin chain tips got the ability to establish bonding with nitrogen atoms of the curative molecular ends that meet a specified distance criteria which is 5 Å. The main reason for choosing this distance is to prevent the carbon atom from establishing a bond with a far distant atom which would lead to numerical difficulties in the energy minimization at the preprocessing stage due to high computational efforts. Using LAMMPS as a solver, 50000 steps of NVE ensemble are applied on the system after the energy minimization. These steps give the chains a chance to cross-link with the curative molecular and distributed homogeneously through the simulation cell. The resulting epoxy network had the cross-link density of 76.3% . Finally, the system is cooled 3.4. Model description and calculations 43 Figure 3.9: The self-avoiding algorithm applied to create the resin and curative chains in space down to room temperature using 50000 NPT steps and relaxed to zero stress using 50000 NPT steps at room temperature. Thereafter, the simulation cell is ready for the numerical experiments figure (3.10). To estimate the thermomechanical properties of the epoxy/clay nano-composites, a new simulation cell is prepared by inserting a (100 × 100 × 10)Å silicate layer surrounded by the organic modifier chains figure (3.11) inside the epoxy box. Furthermore, 50000 steps of NVE are applied on the simulation cell after the energy minimization. These steps induce the epoxy chains to start entanglement 3.4. Model description and calculations 44 Figure 3.10: Simulation cell of the pure epoxy polymer for the molecular dynamics model Figure 3.11: Molecular structure of the silicate layer with addition of the organic modifier’s chains with the om chains that surround the silicate layer. Afterward, the epoxy/clay nano- composite box is cooled down to room temperature and relaxed to zero stress con- figuration using 50000 NPT steps. Up to now, we have two simulation cells. One for the (pure) epoxy polymer and the other for epoxy/clay nano-composite, figure (3.12). These simulation cells are used to estimate the elastic constants as well as the transition temperature and thermal expansion coefficients. The main target of this study is to determine the full parameters of the debonding process at the interface between the epoxy matrix and the silicate layers. Moreover, the orientation of the silicate layer with respect to the external uniaxial load has great effects on the debonding process parameters. Such parameters are hard to measure using the experimental methods. Thus, three simulation cells are prepared not only to estimate the debonding parameters but also to investigate the effects of the silicate 3.4. Model description and calculations 45 Figure 3.12: Simulation cell of epoxy/nano-clay composites for the normal fracture test layer orientations. Starting from the normal debonding where the silicate layer is perpendicular to the external uniaxial load, a (100 × 100 × 250) Å simulation cell has been created. The model contains two parallel silicate layers at each end in the z-direction. The silicate layers were built in the xy-plane and they are separated by 350 chains of the epoxy matrix and organic modifier figure (3.13). An external load was applied in z- direction to investigate the stress-separation curve and the nucleation process around the silicate layer. On the other hand, the shear debonding process is modeled using the simulation cell shown in figure (3.14 left). Here the epoxy matrix must be deformed parallel and opposite to the silicate layer longitudinal direction. Contrasting from the previous simulation cell, the epoxy matrix was divided into two symmetric parts by a gap in x-direction. The width of the inserted gap is selected to be equal to the maximum interaction distance between any two atoms to eliminate any loss of energy and also to concentrate the external load on the interface area. And what is more, the silicate layer lies down in the middle of the simulation cell parallel to the applied loading direction. 3.4. Model description and calculations 46 Figure 3.13: Approximate representation of the nano-clay layers and epoxy polymer in a simulation cell Between the normal and shear debonding, all the parameters of the cohesive model can be approximated by a linear resultant. Nevertheless, an additional test to an orientation in between can show the real behavior of the materials at the interface zone. In order to achieve such a test, a simulation cell similar to that of the shear test has been created. There is only one different, the silicate layer tilts the external load axes by 45◦ as shown in figure (3.14 right). Figure 3.14: Left: Simulation cell of the epoxy/nano-clay composites for the shear mode fracture test. Right: for the test in which the silicate layer is aligned in slope 45◦ with respect to the direction of the external load 3.4. Model description and calculations 47 3.4.3 Calculation of Elastic Constants After achieving the thermomechanical equilibrium of the ensemble of atomistic struc- ture, many methods are available to calculate the elastic constants. Majority liter- ature results present two main methods, namely, stress fluctuation and strain fluc- tuation. In this work, stress fluctuation (static) method, which was developed by Theodorou and Suter [167] is used. Starting from the energy approach, the total potential energy Upot min,0 expression can be extended using a Taylor expansion in the degree of deformation ϵ as: Upot min = Upot min,0 + V0 ∑ ij [τij + ρ0cTγij] ϵij + 1 2 V0 ∑ ij ∑ kl Cijklϵijϵkl (3.61) where the second term in eq. (3.61) represents the internal stress σint such that σint = τ + ρ0cTγ (3.62) where, ρ0 is the density, c is the heat capacity per unit mass, τ is the stress tensor and γ is the Gruneisen tensor [183] which equals to γij = 1 ρ0c αζ kT (3.63) where αζ is the volumetric expansion coefficient. In order to calculate all the elements of the fourth order tensor Cijkl, three loading cases are applied to the undeformed configuration. The first is the uniform hydrostatic compression where ϵ = −ϵ/3 is applied on the diag. (1,1,1) such that ϵ = [−ϵ/3 − ϵ/3 − ϵ/3 0 0 0]T (3.64) Applying such a loading case brings equation (3.61) to Upot min = Upot min,0 − V0 ( −P + αζT kT ) ϵ+ 1 2 V0Bϵ 2 (3.65) 3.4. Model description and calculations 48 where B = C11 + C22 + C33 + 2C23 + 2C31 + 2C12 ∂Upot min ∂ϵ = V0 ( −P + αζT kT ) (3.66) ∂2Upot min ∂ϵ2 = V0B (3.67) Secondly, a pure uniaxial tension along the x-axis is applied as below ϵ = [ϵ 0 0 0 0 0]T (3.68) After submission, equation (3.61) will give the following results for the C11 element of the Cijkl tensor, ∂2Upot min ∂ϵ2 = V0Cii i = (1, 2, 3) (3.69) Similarly, applying the same load again along y and z-axis gives the same formula for C22 and C33 respectively. Finally, a pure shear load is required to estimate the elements C44, C55, and C66. For the zx-shear transformation, the vector of strain components is: ϵ = [0 0 0 0 ϵ 0 ]T (3.70) The second derivative of the Upot min is illustrated as ∂2Upot min ∂ϵ2 = V0Cjj j = (4, 5, 6) (3.71) The main steps of the mentioned method are demonstrated in figure (3.15) below. 3.4.4 Calculation of the interfacial strength parameters The interaction between the internal surface of the epoxy polymer and the external surface of the silicate layer depends on the bilinear cohesive zone model which is based on a traction-separation law. The model is developed by Hillerborg et al. [75]. The 3.4. Model description and calculations 49 Start from undeformed configuration Apply a very small strain to one of the strain tensor components ϵij Minimize the total energy of the deformed configuration Calculate the stiffness matrix component as: Cij = σtens.−σcomp. 2ϵij Calculate the stress as: σ = −Etot V0 Figure 3.15: Stress fluctuation (static) method used to determine the elastic constants traction-separation law, figure (3.16), can be divided into two main zones. The first zone is the linear elastic part 0 ≤ δ ≤ δinit . The constitutive equation is given by σ = Kn ∗ δ (3.72) where Kn is the interphase stiffness. The second zone is the damage initiation and evolution δinit ≤ δ ≤ δfail governed by the following expression σ = σmax ( δfail − δ δfail − δinit ) (3.73) 3.5. Results and discussion 50 σ δ Kn δinit δfail σmax G Figure 3.16: Traction-separation curve at the interfacial zone Therefore, the area below the traction-separation curve represents the fracture toughness G = 1 2 σmaxδfail which is the dissipated energy after complete debonding between the polymer matrix and the silicate layer. 3.5 Results and discussion 3.5.1 Elastic Constants The orthotropic nature of the pure polymer reduces the independent elastic constants of the 4th order stiffness matrix to 9 elements. After applying the constant strain method, the resulting elastic constants for the pure polymer are shown in table (3.1) below C11 C22 C33 C12 C13 C23 C44 C55 C66 3.28 3.20 3.09 1.31 1.26 1.27 0.81 0.79 0.74 Table 3.1: The resulting components of the symmetric 4th order stiffness matrix for the pure epoxy polymer molecular dynamics model The estimated Young modulus, shear modulus, and Poisson ratio are illustrated in table (3.2). Furthermore, a comparison with experimental results has been made. The numerical results have shown fairly good agreement with the experimental mea- surement. 3.5. Results and discussion 51 Young’s Modulus GPa Shear Modulus GPa Poisson’s Ratio Exp. 2.05 0.76 0.35 Sim. 2.52 0.98 0.28 Table 3.2: Mechanical properities of epoxy comparison with experimental data [34] For the epoxy/clay nano-composite, the resulting elastic constants are presented in table (3.3). It is clear that the results show a small variation in z-direction for both C33 and C66. This is due to the high difference in thickness in z-direction relative to the thickness in xy-plane. In other words, the mechanical strength of the simulation cell in z-direction highly depends on the epoxy matrix because of the small thickness of the silicate layer. Similarly, the shear strength of the simulation cell in the xy-plane shows a relatively higher value in comparison with those in the xz- and yz-plane. C11 C22 C33 C12 C13 C23 C44 C55 C66 22.53 22.97 20.50 8.51 7.93 8.06 2.90 3.31 5.09 Table 3.3: The resulting components of the symmetric 4th order stiffness matrix in MPa for the epoxy/nano-clay composites molecular dynamics model At the atomistic scale and for epoxy/clay nano-composite, comparing the numer- ical results with the experimental results is very difficult. Instead of that, we can use elastic constants of a single lamella of MMT (which estimated by Manevitch and Rut- ledge [104]) along with the mixture rule [12] to predict analytical values to compare with. Table (3.4) contains the mechanical properties of epoxy/clay nano-composite in comparison with analytical results. 3.5.2 Transition temperature and thermal expansion coeffi- cient The transition temperatures of both the pure epoxy and epoxy/clay nanocomposite have been tested by applying linear increments of temperature. The simulation cell develops a linear growth in its volume with respect to the increment of the applied 3.5. Results and discussion 52 (GPa) Young’s Mod. (GPa) Shear Modulus Poisson’s Ratio Ex Ey Ez Gyz Gxz Gxy Analytical 21.944 21.947 21.441 3.453 3.457 7.676 0.41 Simulation 17.865 17.864 16.393 2.891 3.283 5.071 0.27 Table 3.4: Mechanical properties of epoxy/nano-clay composites in comparison with analytical data temperature. During the simulation, the volume-temperature curve shows an abrupt slope at a specific region. The resulting curve shows two clear linear behavior zones. The transition temperature Tg can be estimated by intersecting the two linear zones as shown in figure (3.17). The estimated Tg of the pure epoxy matrix is 375 K which is still inside the experimental range. Furthermore, the thermal expansion coefficients of both the glassy and the rubbery phase are calculated from the slope of the volume- temperature curve. Temperature K Vo lu m e n m 3 Figure 3.17: Volume-Temperature curve for the pure epoxy matrix Due to the huge variation between the thermal expansion coefficient of the epoxy matrix and the silicate layer, the epoxy/clay nano-composite simulation cell expresses 3.5. Results and discussion 53 a higher transition temperature (Tg=393 K) in comparison with the pure epoxy cell as illustrated in fig (3.18) below. Moreover, insertion of a silicate layer inside the simulation cell increases the crosslinking density by reason of the chains entanglement between the epoxy from one side and the organic modifier from another side. As it is known, the crosslinking density has linear direct proportion with Tg as presented by [92, 151, 158]. Vo lu m e n m 3 Temperature K Figure 3.18: Volume-Temperature curve for the epoxy/nano-clay composites 3.5.3 Fracture characterization of the epoxy/clay interface As reported by many literatures, the strength of debonding and void nucleation in the interphase zone between the matrix and inclusions extremely affects the ther- momechanical properties of the material composite at the micro and macro levels. Furthermore, the void nucleation and debonding processes are initiated early during the elastoplastic behavior under external loads. In this section, the Molecular Dynamics (MD) method is used to estimate stress- separation curves between the epoxy matrix and the silicate layers. In addition, the effects of the silicate layer orientations are tested for three different orientations. The 3.5. Results and discussion 54 first orientation is the normal, in which the silicate layer set in a position perpendic- ular to the external uniaxial load direction figure (3.19). Epoxy polymer Silicate layer Figure 3.19: Molecular configuration of epoxy/nano-clay interface during normal debonding simulation (external load is perpendicular of the silicate layer) The external load is resisted by two main internal forces. The first one is the cohesion of the polymer matrix which results from the cross-link density and the entanglement between the polymer chains. The second one is the adhesion which can be considered as a summation of the interactions at the main two interfacial areas: the interactions between the nano-clay layer and the organic modifiers (due to the strong ionic bond that is created during the exchange reaction of the modification process [18]) and the interaction between the epoxy polymer and the organic modifiers (this interaction results from the weak van der Waals interactions in addition to the entanglements between the chains of the two materials). As a result of the high cross- link density of the epoxy polymer, the cohesive force is found to be much higher than the adhesive force. This can explain why debonding takes place at the interfacial area between the epoxy polymer and the organic modifier. Figure (3.20) shows the traction-separation curve for the epoxy polymers filled with nano-clay layers. Three main parameters can be estimated from the resulting curve. Namely, the fracture toughness, damage initiation separation, and damage failure separation. They are necessary for the next scale (Microscale) of the model. Moreover, another parameter can be investigated from the curve. The fracture energy, which can 3.5. Results and discussion 55 be calculated by finding the area under the traction-separation curve, is considered as the critical energy at which the crack propagation is going to start. St re ss (M Pa ) Seperation (nm) Figure 3.20: Stress-separation curve during normal debonding at 300 K The sudden disentanglement between the matrix and organic modifier chains ap- pear as fast drop in stress in the tension separation curve after the damage initiation. That explains the big role that the chains entanglement between epoxy and the or- ganic modifier plays during the debonding process. In contrast, when the silicate layer is parallel to the external uniaxial load direction (shear debonding), the debonding process is illustrated in figure (3.21) below. Fur- thermore, the stress-separation curve, figure (3.22) shows relatively lower parameters in stress and higher values in separation. Between normal and shear debonding, a linear variation of the cohesive parame- ters with respect to the silicate layer orientation can be assumed. Nevertheless, an additional numerical test has been performed with silicate layer sets at an angle of 45◦ with respect to the external uniaxial load direction. The results are shown in figure (3.23) the corresponding tension separation curve is presented in figure (3.