Michigan Technological University

Digital Commons @ Michigan

Tech
Dissertations, Master's Theses and Master's Reports
Dissertations, Master's Theses and Master's Reports
- Open

2015

MOLECULAR MODELING OF EPON

862-DETDA / CARBON COMPOSITES
Cameron Hadden
Michigan Technological University

Copyright 2015 Cameron Hadden

Recommended Citation
Hadden, Cameron, "MOLECULAR MODELING OF EPON 862-DETDA / CARBON COMPOSITES", Dissertation, Michigan
Technological University, 2015.
http://digitalcommons.mtu.edu/etds/946

Follow this and additional works at: http://digitalcommons.mtu.edu/etds

Part of the Mechanical Engineering Commons

MOLECULAR MODELING OF EPON 862-DETDA / CARBON
COMPOSITES

By
Cameron Hadden

A DISSERTATION
Submitted in partial fulfillment of the requirements for the degree of

DOCTOR OF PHILOSOPHY
In Mechanical Engineering – Engineering Mechanics

MICHIGAN TECHNOLOGICAL UNIVERSITY

2015

© 2015 Cameron M. Hadden

This dissertation is hereby approved in partial fulfillment of the requirements for
the Degree of DOCTOR OF PHILOSOPHY LQ MHFKDQLFDO EQJLQHHULQJ-
EQJLQHHULQJ MHFKDQLFV.

Department of Mechanical Engineering – Engineering Mechanics

Dissertation Advisor: Dr. Gregory M. Odegard

Committee Member: Dr. Reza Shahbazian Yassar

Committee Member: Dr. Julia King

Committee Member: Dr. Gowtham

Department Chair: Dr. William Predabon

 Table of Contents

List of Figures .......................................................................................................... vii

List of Tables ............................................................................................................... x

Preface ...................................................................................................................... xii
Acknowledgments ................................................................................................... xiii
Abstract ..................................................................................................................... xv

Chapter 1: Introduction ............................................................................................. 1

1.1 Material Overview .......................................................................................... 1

Epoxy Resins......................................................................................................... 1

Carbon Reinforcements...........................................................................................

1.2 Molecular Dynamics ....................................................................................... 4

Chapter 2: Molecular Modeling Procedures ........................................................... 7

2.1 Force Field ....................................................................................................... 7

2.2 Modeling Initial Structures .......................................................................... 9

2.3 Crosslinking Procedure ............................................................................... 12

Chapter 3: Modeling the Fiber/Matrix Interface .................................................. 16

3.1 Introduction .................................................................................................. 16

3.2 Modeling Procedure .................................................................................... 17

3.3 Results and Discussion ................................................................................ 22

v
 Polymer Mass Density Characteristics ............................................................... 23

Polymer Internal Stress Characteristics .............................................................. 26

Polymer Potential Energy Characteristics .......................................................... 28

3.4 Conclusions .................................................................................................... 29

Chapter 4: Multi-Scale Modeling of Graphene Doped Epoxy ............................. 30

4.1 Introduction .................................................................................................. 30

4.2 Modeling Procedure ..................................................................................... 31

Mechanical Deformation .................................................................................... 35

Micromechanics.................................................................................................. 37

4.3 Results and Discussion ................................................................................. 39

4.4 Conclusions .................................................................................................. 46

Chapter 5: Multi-Scale Modeling of Graphene Oxide / Epoxy Composites ....... 47

5.1 Introduction .................................................................................................. 47

5.2 Modeling Procedure ..................................................................................... 49

5.3 Results and Discussion ................................................................................. 53

5.4 Conclusions .................................................................................................. 57

Chapter 6: Scope for Future Works ....................................................................... 58

6.1 Improving Modeling Procedures ................................................................ 58

6.2 Modeling Functionalized Graphene and Surface Treatments ................. 59

6.3 Modeling Aging in Composite Structures .................................................. 60

vi
References ................................................................................................................. 62

Appendix A
LAMMPS scripts for MD simulations, MM simulations and Deformation
Simulations ................................................................................................................ 68
A.1 Molecular Energy Minimization script ............................................................ 68

A.2 MD NVT simulation script for equilibration ................................................... 69

A.3 MD NPT simulation script for tensile deformation ......................................... 72

A.4 MD NVT simulation script for shear deformation........................................... 76

A.5 MD OPLS UA parameters from coordinates file............................................. 77

Appendix B
Copyright Agreements and Permissions ................................................................ 82
B.1 Permission to use Figure 1.2 ............................................................................ 82

B.2 Copyright agreement for using selected text and Figures in Chapter 2 ........... 83

vii
 List of Figures

Figure 1.1 – Epoxed (a) and Amine (b) moleculuar structures .................................... 2

Figure 1.2 – Predicted internal structure for carbon fibers, surface atom structures,
and aromatic carbon ring structure ............................................................................... 3

Figure 2.1 – A visual representation of each energy contribution that is considered by

the OPLS United Atom force field............................................................................... 8

Figure 2.2 – Molecular structures of a) EPON 862 and b) DETDA shown with
simulated image ......................................................................................................... 10

Figure 2.3 – Visual representation of the fully equilibrated 2:1 ratio of EPON
862:DETDA ............................................................................................................... 12

Figure 2.4 – Diagram of crosslinking reaction .......................................................... 13

Figure 2.5 – Flowchart describing crosslinking algorithm ........................................ 14

Figure 3.1 – Resulting structure after implementing the ‘randomly rotate and
translate’ program for multiplying the initial 2:1 ratio structure ............................... 18

Figure 3.2 – a) Initial planar graphite structure and b) Graphite structure after
equilibration ............................................................................................................... 19

Figure 3.3 – Combining the three separate MD models to form the complete MD
model of the composite interface ............................................................................... 21

Figure 3.4 – Crosslink atom mass density along the z-axis for each crosslinked
structure ...................................................................................................................... 22

Figure 3.5 – Mass density profile along Z-axis of simulation box for graphite and
polymer molecules ..................................................................................................... 23

Figure 3.6 – Zoomed in interfacial mass density profile along Z-axis ...................... 24

Figure 3.7 – TEM image of carbon nanofiber embedded in an epoxy resin ............. 25

viii
Figure 3.8 – EPON 862 and DETDA molecular mass densities along simulation box
Z-axis for Non-crosslinked and 80% crosslinked structures ...................................... 26

Figure 3.9 – Spatially averaged Von Mises stress for polymer and graphite atoms
along Z-axis ................................................................................................................ 27

Figure 3.10 – Spatially averaged per-atom potential energy along Z-axis................. 28

Figure 4.1 – Molecular structures for single layer graphene sheet model highlighting
the densification process............................................................................................. 33

Figure 4.2 – Equilibrated models for varying number of graphene layers before
crosslinking................................................................................................................. 34

Figure 4.3 – Crosslink dispersion as a function of z-axis coordinate for 80%

crosslinked systems .................................................................................................... 35

Figure 4.4 – Simulated deformation characteristics for tension along the z-axis and
shear in the x-z direction ............................................................................................ 36

Figure 4.5 – Tensile and Shear stress vs strain curves for varying number of
graphene layers in 80% crosslinked structures ........................................................... 36

Figure 4.6 – Schematic representation for implementing MD results into MACGMC

software to obtain data for GNP doped epoxy and GNP epoxy/Carbon Fiber
composites .................................................................................................................. 38

Figure 4.7 – MAC/GMC RUC for GNP doped epoxy / carbon fiber model ............. 39

Figure 4.8 – Atom density profile along z-axis for a)1 layer b) 2 layers c) 3 layers
and d) 4 layers of graphene ....................................................................................... 40

Figure 4.9 – Atomic stress before and after tensile deformation along the y-axis. .... 41

Figure 4.10 – MAC/GMC data for addition of GNPs with carrying number of layers
and interfacial polymer crosslink density ................................................................... 43

Figure 4.11 – MAC/GMC tensile modulus for GNP doped epoxy/carbon fiber
models for varying GNP volume fractions and atomic thicknesses .......................... 44

ix
Figure 4.12 – Experimental and MAC/GMC tensile modulus for GNP doped
epoxy/carbon fiber models ......................................................................................... 45

Figure 4.13 – MAC/GMC transverse modulus for GNP doped epoxy/carbon fiber
models for varying GNP volume fractions and atomic thicknesses .......................... 46

Figure 5.1 – Widely accepted models for the structure of graphene oxide................ 48

Figure 5.2 – The most stable confirmations of graphene oxide as found by DFT for
a) epoxide oxygen only, b) hydroxyl molecules only, and c) a combination of
epoxide and hydroxyl molecules ................................................................................ 51

Figure 5.3 – Single and 4 layer graphene oxide structures containing 12% bonded
oxygen species ........................................................................................................... 51

Figure 5.4 – Fully equilibrated Single and 4 layer graphene oxide / epoxy structures
containing 12% bonded oxygen species .................................................................... 52

Figure 5.5 – Polymer mass density profile viewed along the Z-Axis for single layer
a) 12% bonded oxygen b) 24% bonded oxygen and 4 layer c) 12% bonded oxygen d)
24% bonded oxygen .................................................................................................. 54

Figure 5.6 – Visual representation of per-atom potential energy values (Kcal/mol-

Angstrom) for single and 4 layer graphite oxide structures for varying oxygen
species with 12% oxygen bonding ............................................................................ 55

x
 List of Tables

Table 1.1 – Common carbon fibers and their corresponding mechanical properties.. 3

Table 4.1 – MAC/GMC Input Parameters for HexTow® AS4-GP/3K Fibers .......... 39

Table 4.2 – Resulting data from MD simulations....................................................... 41

Table 4.3 – Experimental results for epoxy/carbon fiber composites with varying
GNP volume fractions ................................................................................................ 45

Table 5.1 – Mechanical properties derived from MD for oxidized graphene models
along with those for 75% crosslinked pristine graphene MD models from chapter 4
for comparison ........................................................................................................... 56

Table 5.2 – MAC/GMC results for 75% corsslinked pristine and epoxide oxidized
graphite structures of 1 and 4 atomic layers ............................................................... 57

xi
 PREFACE

The work outlined in this dissertation was carried out in the department of
Mechanical Engineering – Engineering Mechanics at Michigan Technological
University over the period of August 2009 to December 2014. Portions of Chapters 2
and 3 have been presented in the following publication:
C.M. Hadden, B.D. Jensen, A. Bandyopadhyay, G.M.Odegard, A.Koo, R.Liang
Composite Science and Technology 2013, 76, 92-99. © Elsevier Ltd
Portions of chapter 4 have been submitted for publication, and the information
outlined in Chapter 5 could potentially be used for publication in the near future.
All of the writing in this dissertation was performed by yours truly, as well as all of
the Molecular Dynamics research, which includes creating and deforming the models
for extraction of material properties.
For help from other parties, I would first like to thank Dr. Gregory Odegard for his
help with writing and formatting journal article submissions. I would like to thank
Ananyo Bandyopadhyay for his molecular modelling contributions in Chapter 3, as
well as Ben Jensen for help with programming. Also from Chapter 3, I would like to
thank A. Koo and R. Liang for the use of their TEM image, which helped to
reinforce my research. I would like to thank Dr. Julie King and Danielle Klimek for
their contributions to Chapter 4, which include fabrication and mechanical testing of
carbon composites for comparison to simulation results. Finally, I would like to
thank Dr. Evan Pineda from NASA research for developing specific modelling
scripts that would enable simulation data from Chapters 4 and 5 to be input into
micromechanics software.

xii
 ACKNOWLEDGMENTS

Obtaining my Ph.D. has been a great experience. I have been able to learn a lot about
mechanical engineering, specifically mechanics of composite materials. I have also
learned to do some pretty intricate programming and taken a few courses in materials
science. I am very thankful for having such a great advisor in Dr. Gregory Odegard.
He has helped to steer me through my doctorate and put me in a great position to
obtain the job I wanted when accepting the challenge of obtaining my Ph.D.
I would also like to thank Dr. Reza Shahbazian-Yassar, Dr. Tolou Shokuhfar, Dr.
Ranjit Pati, Dr. Gowtham Shankara, and Dr. Julia King for agreeing to be on my
Ph.D. committee. The members of my committee changed from my initial proposal
defense to my final defense, though all the aforementioned names have supported my
research and provided valuable insight. During my Ph.D. I was able to travel and
work with several research groups. From Dayton Air Force base I would like to
thank Dr. Ajit Roy and Dr. Vikas Varshney. From NASA Langley I would like to
thank Dr. Kris Wise. From NASA Glenn I would like to thank Dr. Evan Pineda and
Dr. Brett Bednarcyk. All of these gentleman have put in a lot of hard work in relation
to Molecular Dynamics and micromechanics, without their guidance things would
have been much more difficult for me.
This research would not have been possible without the funding and support from
both NASA and the Air Force. Funding was provided by NASA under the Aircraft
Aging and Durability Project (Grant NNX07AU58A) and the Revolutionary
Technology Challenges Program (Grant NNX09AM50A); and by the Air Force
Office of Scientific Research under the Low Density Materials Program (Grant
FA9550-09-1-0375).
From Michigan Technologicaly University, I would like to thank the other members
of my research group. Ananyo Bandyopadhyay helped me tremendously in getting a
grasp on Molecular Dynamics and setting up my initial systems. Ben Jensen
provided me with a lot of help with the various programs I had to write to create my
models. I’d also like to thank Dr. Gowtham Shankara for all of his help with the
cluster and his programming classes. All of those previously mentioned have been
both great co-workers and great friends.
This has truly been an amazing experience for me. Houghton has not only provided
me with a superior education, but has enabled me to enjoy the extra-curricular
activities that make me who I am. I have been blessed to snowboard in the deepest
powder, wakeboard on the flattest water, and surf some amazing waves on Lake
Superior. I’d like to thank all of those outside of the university who were there to
share in the excitement and for all of the high fives. I have been able to maintain the
perfect balance of work and play, and I’ve had great people by my side for all of it.
xiii
Finally, I would like to thank my family and friends that have been there from the
beginning. My parents were the driving force behind me choosing Michigan Tech.
They supported me throughout all the decisions I made while I was here, including
that first summer spent away from home. My brother and sister have always
encouraged me to do my best, and have been quick to tell others how proud they are
of me. I will never forget that. Last but not least, I would like to say thank you to my
lovely girlfriend, Jaime LeDuc, who has been there through all the ups and downs.
She has been my most powerful inspiration for success.

Cameron Hadden
August 2014
Houghton, MI
USA

xiv
 ABSTRACT
The thermoset epoxy resin EPON 862, coupled with the DETDA hardening agent,
are utilized as the polymer matrix component in many graphite (carbon fiber)
composites. Because it is difficult to experimentally characterize the interfacial
region, computational molecular modeling is a necessary tool for understanding the
influence of the interfacial molecular structure on bulk-level material properties. The
purpose of this research is to investigate the many possible variables that may
influence the interfacial structure and the effect they will have on the mechanical
behavior of the bulk level composite. Molecular models are established for EPON
862-DETDA polymer in the presence of a graphite surface. Material characteristics
such as polymer mass-density, residual stresses, and molecular potential energy are
investigated near the polymer/fiber interface. Because the exact degree of cross-
linking in these thermoset systems is not known, many different crosslink densities
(degrees of curing) are investigated. It is determined that a region exists near the
carbon fiber surface in which the polymer mass density is different than that of the
bulk mass density. These surface effects extend ~10 Å into the polymer from the
center of the outermost graphite layer. Early simulations predict polymer residual
stress levels to be higher near the graphite surface. It is also seen that the molecular
potential energy in polymer atoms decreases with increasing crosslink density.
New models are then established in order to investigate the interface between EPON
862-DETDA polymer and graphene nanoplatelets (GNPs) of various atomic
thicknesses. Mechanical properties are extracted from the models using Molecular
Dynamics techniques. These properties are then implemented into micromechanics
software that utilizes the generalized method of cells to create representations of
macro-scale composites. Micromechanics models are created representing GNP
doped epoxy with varying number of graphene layers and interfacial polymer cross-
link densities. The initial micromechanics results for the GNP doped epoxy are then
taken to represent the matrix component and are re-run through the micromechanics
software with the addition of a carbon fiber to simulate a GNP doped epoxy/carbon
fiber composite. Micromechanics results agree well with experimental data, and
indicate GNPs of 1 to 2 atomic layers to be highly favorable.
The effect of oxygen bonded to the surface of the GNPs is lastly investigated.
Molecular Models are created for systems with varying graphene atomic thickness,
along with different amounts of oxygen species attached to them. Models are created
for graphene containing hydroxyl groups only, epoxide groups only, and a
combination of epoxide and hydroxyl groups. Results show models of oxidized
graphene to decrease in both tensile and shear modulus. Attaching only epoxide
groups gives the best results for mechanical properties, though pristine graphene is
still favored

xv
 CHAPTER 1: INTRODUCTION

1.1 Material Overview

Thermoset epoxy/carbon fiber composites have long been one of the primary
material choices for modern aerospace applications. These materials exhibit an
excellent strength to weight ratio and can be easily manufactured to have specific
performance characteristics. The overall material properties of epoxy composites
depend not only on the individual properties of both the carbon fiber and epoxy
matrix, but also the physical and chemical conditions of the fiber/matrix interface.1 A
good bond between the fiber and matrix enables improved load transfer from the
matrix to the fibers which results in greater composite stiffness and strength, but
lower toughness because of the absence of crack deflection mechanisms at the
interface. In contrast, a weak bond between fiber and matrix will provide lower
composite strength with a higher toughness.2 The molecular behavior of the epoxy
matrix and carbon fiber surface can provide valuable insight to the interface, and
thereby the mechanical behavior of the entire composite.

Epoxy Resins

Thermoset epoxy resins derive their name from the 3-membered epoxide ring
found in their molecular structure, as shown in Figure 1.1(a). The epoxide ring is a
high energy molecular structure and this energy is reduced by reacting to form bonds
with different molecules, commonly known as hardeners. Initially both the epoxy
and the hardener components are a viscous liquid, and when mixed they polymerize
and harden. This process is often referred to as curing. It is important to note that the
extent of curing in thermoset epoxies is unknown, and will be an important variable
for the research performed herein. Different epoxy resins can have different amounts
of epoxide groups per molecule, and different epoxy molecules will produce a
polymer with unique mechanical properties. Hardeners often contain a reactive

1
hydrogen group such as an amine, shown in Figure 1.1(b), or acid. Different
hardeners will also produce a polymer with unique properties.3

Figure 1.1 a) Epoxide ring containing two reactive carbon atoms and one oxygen atom. b) Amine
group containing two reactive hydrogen atoms and one nitrogen atom.

Fully cured epoxy resins have many desirable properties such as: excellent
adhesion to many materials, good mechanical and electrical properties, high heat
resistance, and resistance to chemicals and moisture.3 Because of this they have
many applications, from aerospace composites to house-hold adhesives.

Carbon Reinforcements

The most common form of carbon reinforcement used in modern composites

are carbon fibers. These fibers are usually manufactured by the
oxidation/carbonization of polyacrylonitrile (PAN) fibers.1 Carbon fibers can be
manufactured to exhibit an array of material properties, and are commonly
categorized as either high modulus (Type I) or high strength fibers (Type II).4 High
modulus fibers UHTXLUHKHDWWUHDWPHQWGXULQJPDQXIDFWXULQJRI!ࡈ&ZKLOHKLJK
VWUHQJWKILEHUVDUHPDQXIDFWXUHGDWaࡈ&4 The mechanical properties for carbon
fibers depend strongly on the manufacturing process and the test method used to
determine the desired property. A list of common carbon fibers and there
corresponding mechanical properties is provided in Table 1.1. Carbon fibers are used
in composites by weaving individual fibers into a desirable pattern and combining
them with a polymer matrix. Figure 1.2 depicts a highly accepted structure of a
carbon fiber. It can be seen that the surface molecules are much more organized that
that of the core, and the molecules on the surface are overlapping sheets composed of
aromatic carbon rings.
2
Table 1.1 Common carbon fibers and their corresponding mechanical properties4

Figure 1.2 Predicted internal structure for carbon fibers, surface atom structures, and aromatic carbon
ring structure. Image courtesy of NASA Report 4084 © L.T. Drzal

Surface treatments are often applied to carbon fibers to improve bonding to

matrix polymers.5 Chemical treatments can increase the surface energy of carbon
fibers and promote bonding to matrix polymers by adding reactive species that are

3
able to bond to the certain groups within the matrix polymer molecules. Fibers can
also be coated with ‘sizings’ to protect them during the fabrication process, and in
some cases increase bonding.11 Carbon nanostructures such as nanotubes6,7 and
graphene nanoplatelets (GNPs)8,9 can also be used within a polymer matrix to create
nano-composites. They have also been shown to increase mechanical properties
when added to carbon fiber/epoxy composites.7,10

1.2 Molecular Dynamics

Molecular Dynamics (MD) can be defined as a form of computer simulation

in which atoms and molecules interact for a defined period of time by
approximations of known physics, thus describing the motion of the atoms. MD is
based on statistical mechanics, and requires a lot of computational power. There are
four types of MD simulations most commonly used:12

1. Canonical ensemble (NVT ensemble)- The NVT ensemble does not allow the
system to undergo changes in number of moles (N), volume (V), or
temperature (T). A thermostat is needed to help steer the exchange of
exothermic and endothermic processes within the system. Different
thermostats are available to remove energy from the system in a manner
‘realistic’ for specific applications.

2. Isothermal-Isobaric (NPT ensemble)- The NPT ensemble does not allow for
change to number of moles (N), pressure (P), or temperature (T). A
thermostat and barostat are necessary to help control temperature and
pressure, respectively. The volume is allowed to fluctuate based on the
temperature and pressure of the system.

4
 3. Microcanonical ensemble (NVE ensemble)- The NVE ensemble does not
allow the system to undergo changes in number of moles (N), volume (V), or
energy (E). By keeping the energy term constant, it creates an adiabatic
process and does not allow for heat exchange. Overall energy is kept constant
and is continuously changing between kinetic and potential energies.

4. Isoenthalpic-Isobaric (NPH ensemble)- The NPH ensemble does not allow

change to number of moles (N), pressure (P), or enthalpy (H). This ensemble
denotes volume (V) as a dynamic variable relating to the coordinates of all
the atoms in the system. Kinetic energy is defined as PV. Enthalpy can then
be defined as H = E + PV, where H is kept constant and internal energy (E)
and kinetic energy PV are allowed to change in response to one another.

Energy minimization, or molecular minimization (MM), simulations are used

in conjunction with MD to reduce the overall energy of a system. MM simulations
make small adjustments to molecular geometries searching for low energy
confirmations. The ultimate goal in equilibrating an MD model is to search for the
lowest possible energy confirmation, as that is what occurs in nature based on the
laws of thermodynamics.

The use of computational methods for studying material behavior has really
taken off in the last decade. MD gives researchers insight to the molecular behavior
of a material, something that experimental techniques have yet to explain. The most
disputed aspect of MD is the timescale for which the simulations are conducted.
Most MD simulations occur over a few hundred picoseconds, while longer
simulations may span a few nanoseconds. Because of this, they are difficult for
comparison with experimental data. As MD modeling techniques continue to
improve, they will enable researchers to develop new materials outside of the lab
setting and will provide information that will steer future experimental research. By
teaming up with experimental researchers, MD can help lessen time in developing
new materials and the money associated with experimental techniques.
5
 Much research has been done using MD for thermoset epoxies. Early MD
methods for studying thermoset epoxies utilized bead-spring models13,14 and Monte-
Carlo simulations based on bond fluctuation models15-17. More recent techniques for
studying epoxies utilize detailed molecular structures and produce models of much
larger size (more atoms). Recent models have also focused heavily on the
crosslinking mechanisms that occur during curing.18-21 These studies have lead the
way to developing models for many epoxy/carbon composite systems. There has
been extensive research performed using MD models to investigate thermoset
polymers containing nanotubes22-28, nanoparticles29-31, and in the presence of a
graphite surface (representing either graphene nanoparticles or the surface of a
carbon fiber).33-37 The interface between epoxy and carbon reinforcement is
investigated in nearly all of the aforementioned references. Mechanical properties
have been derived from MD models for a variety of epoxy/carbon systems.36-40
Although these studies have given valuable information regarding the physical
nature of the interfacial region, including an estimate of material properties on the
molecular level, there has been little effort to implement this molecular data into a
larger scale model that will enable comparison to experimental results for similar
systems.

The objective of this research was to first establish an efficient and accurate
procedure for developing large scale models of EPON 862-DETDA epoxy in the
presence of a carbon fiber surface. The interfacial characteristics of the polymer were
analyzed for multiple crosslink densities (degrees of curing). The initial models were
validated with experimental TEM images for similar systems. The methodology used
to establish the initial model would then be used to create new models of EPON 862-
DETDA epoxy with embedded graphene nanoplatelets (GNPs). Mechanical
properties were extracted using MD techniques and the effects of GNP atomic
thickness and polymer crosslink density were investigated. The mechanical

6
properties derived from the MD models were then implemented into a
micromechanics model based on the generalized method of cells. This
micromechanics model would enable comparison to experimental results for both a
GNP doped epoxy specimen, as well as a GNP doped epoxy / carbon fiber
composite. Lastly, the effect of adding oxygen functional groups to the surface of
graphene nanoplatelets was investigated. MD models were created with GNPs of
multiple atomic thicknesses with a variety of oxygen species attached in multiple
concentrations. For this final study, polymer crosslink density was kept at a constant.
Again, micromechanics models were used for comparison to experimental results.

CHAPTER 2: MOLECULAR MODELING PROCEDURES1

This section describes the procedures used for establishing molecular models.
The Force Field chosen for this research is highlighted, followed by the procedure
for building the initial model of a 2:1 stoichiometric ratio of EPON 862:DETDA that
would be used as the building block for all larger models. The crosslinking
procedure used throughout this research will lastly be discussed.

2.1 Force Field and Simulation Parameters

The force field used in this research was the OPLS United Atom force field
developed by Jorgensen and co-workers.41-44 This force field calculates the total
energy of the system as the sum of all the individual energies associated bond, angle,
dihedral and 12-6 Lennard-Jones interactions. A visual representation for each
energy contribution can be seen in Figure
1
Portions of text in this chapter have been reprinted with permission from C.M. Hadden, B.D. Jensen,
A. Bandyopadhyay, G.M.Odegard, A.Koo, R.Liang Composite Science and Technology 2013, 76, 92-
99. © Elsevier Ltd. Please refer Appendix B for copyright information.
7
Figure 2.1 A visual representation of each energy contribution that is considered by the OPLS United
Atom force field.

The equations for each energy contribution are given as follows. The bond energy
contribution is a simple quadratic equation given as:

‫ܧ‬௕௢௡ௗ௦ = σ௕௢௡ௗ௦ ‫ܭ‬௥ (‫ ݎ‬െ ‫ݎ‬଴ )ଶ (2.1)

Kr is a force constant, r is the distance between the two bonded atoms, and r0 is the
equilibrium bond distance for the two bonded atoms. The energy associated with
bond-angle bending is also a quadratic equation and is given as:

‫ܧ‬௔௡௚௟௘௦ = σ௔௡௚௟௘௦ ‫ܭ‬ఏ (ߠ െ ߠ଴ )ଶ (2.2)

Kύ is a force constant, ύ is the angle between the bonds, and ύ0 is the equilibrium
bond angle. The dihedral energy is calculated as:

௏భ ௏మ ௏య ௏ర
‫ܧ‬ௗ௜௛௘ௗ௥௔௟௦ = [1 + cos ‫ ]׎‬+ [1 െ cos(ʹ‫ ])׎‬+ [1 + cos ͵‫ ]׎‬+ [1 െ
ଶ ଶ ଶ ଶ

cos Ͷ‫]׎‬ (2.3)

where V1, V2, V3 and V4 are coefficients of a Fourier series42,44, ĭ is the dihedral
angle. The van der Waal’s interactions are represented by the 12-6 Lennard Jones
equation:

ଵଶ ଺
ఙ ఙ
‫ܧ‬௜௝ = 4ߝ௜௝ [൬ ೔ೕ ൰ െ ൬ ௥ ೔ೕ ൰ ] (2.4)
௥ ೔ೕ ೔ೕ

Where ı is the equilibrium spacing parameter, and was considered to be the

arithmetic mean of the individual parameters of atoms i and j. İis the well depth
parameter and was considered to be the geometric mean of the values given for
atoms i and j. Here, r denotes the distance between the non-bonded atoms i and j.
8
Charges were not considered in this work, as they require a lot of computational
power to calculate. Also graphite structures are known to be neutral of charge and
will be the main influence on interfacial atomic structure, thereby deeming charges
unnecessary for this particular research.

2.2 Modeling Initial Structures

This section describes the procedure for establishing the 2:1 stiochiometric
model of EPON 862 monomer: DETDA hardener. The procedure for multiplying
this initial small structure into a much larger system is different depending on the
desired simulation data, and hence will be described in the appropriate chapters. The
initial 2:1 ratio was produced in a similar matter for all models, utilizing the
LAMMPS (Large Scale Atomic/Molecular Massively Parallel Simulator)12,45
software package for all of the MM and MD simulations described in the following.

The initial uncrosslinked polymer molecular structure was established using a

procedure similar to that of Bandyopadhyay et al.20, consisting of the EPON 862 (Di-
glycidyl ether of Bisphenol-F) monomer and the DETDA (Diethylene Toluene
Diamine) hardener shown in Figure 2.2. A stoichiometric ratio of 2 molecules of
EPON 862 and 1 molecule of DETDA was placed in a 10 x 10 x 10 Å MD
simulation box with periodic boundary conditions. The 2:1 stoichiometric ratio was
chosen due to the fact that if all possible crosslinks were formed, the polymer would
achieve 100% crosslink density (full curing) allowing the polymer to exhibit the best
possible material properties. The initial atomic coordinates file was written in the
native LAMMPS format with the OPLS United Atom force field, as previously
described. The non-bonded van der Waal’s interactions were modeled with in
interaction cutoff radius of 10 Å. The use of the OPLS United Atom force field
allows for modeling of CH3, CH2, CH, and alkyl groups as single united atoms with
their corresponding masses. The described polymer model utilizes united atoms

9
structures for all applicable groups, except for the C and H atoms in the phenyl rings
for both monomer and hardener molecules along with one CH3 group directly
connected to the phenyl ring of the DETDA molecule. Thus, the use of united atoms
reduced the 2:1 structure from 117 atoms to 83 atoms. The location of each united
atom is show in Figure 2.2, with 31 total atoms in the molecule of EPON 862 and 21
in the molecule of DETDA.

Figure 2.2 Molecular structures of a) EPON 862 and b) DETDA shown with simulated image. Top
chemical structures show united atoms to be circled, corresponding to the green beads in the simulated
image.

Creating a coordinates file requires defining specific atom types, which were
easily located in the OPLS United Atom parameters manual. Bonds were defined
according to the atom types for each bonded atom, again these were easily obtained.
All OPLS United Atom force field parameters used in this work are given in Appendix
A.5. Angle and Dihedral parameters are minimally defined in literature, and therefor
required several assumptions to be made.

Assumptions for angle types:

1. C-C-N was given the same parameters as C-C-C

2. CH2-CH-CH2, CH3-CH-CH2, CH3-CH-CH3, CH-CH-CH3, and CH-CH-CH
all utilized the same parameters

10
 3. CH3-CH2-C was given the same parameters as CH-CH2-CH
4. CH-CH2-O was given the same parameters as CH2-CH-O
5. CH2-O-CH, C-O-CH2, CH2-O-CH3, CH2-O-CH2, and CH-O-CH all
utilized the same parameter set

Assumptions for dihedral types:

1. Ca-Ca-CH2-Ca and CH3-CH2- Ca- Ca were given same parameters as Ca- Ca-
Ca- Ca, where Ca = aromatic carbon
2. All dihedral N atoms in DETDA were considered as C atoms
3. C-O-CH2-CH, CH-O-CH-CH, and CH2-O-CH-CH2 were given same
parameters
4. CH2-CH-CH2-O, C-CH-CH2-O and CH2-CH-CH2-O were given same
parameters

The main contribution to the potential energy of the system is caused by bond
energy, as more energy is required to alter bond distances than to distort an angle or
dihedral angle from equilibrium. In saying this, the assumptions made should not
significantly affect the overall potential energy of the model.

After establishing an atomic coordinates file, the initial 2:1 structure was
subjected to four MM minimizations and three MD simulations in order to minimize
internal forces and reduce internal residual stresses created by the initial construction
of bonds, bond angles, and bond dihedrals. The MM simulations were run with a
conjugate gradient stopping criterion which utilizes stored information from
successive line searches to steer the atomic configuration to a minimal energy
confirmation.12 At each iteration the force gradient is combined with information
from the previous iteration to compute a new search direction perpendicular
(conjugate) to the previous search direction. The minimization ceases when no
moves can be made to further reduce the energy of the system. The MD simulations
utilize the NVT ensHPEOHZLWKHDFKVLPXODWLRQODVWLQJSVDWࡈ. A timestep of
0.2 femtoseconds was utilized for initial equilibration steps and all other MD
simulations described in this work. A Nose/Hoover thermostat and barostat was
11
implemented for temperature and pressure control, respectively.46 After equilibrating
the initial 2:1 structure, it was ready to be inserted into a program that would
randomly rotate and translate it to fill a desired volume. The equilibrated 2:1
structure can be seen in Figure 2.3

Figure 2.3 A visual representation of the fully equilibrated 2:1 ratio of EPON 862:DETDA.

2.3 Crosslinking Procedure

The equilibrated model was crosslinked based on the root mean square
(RMS) distance between the CH2 groups of the EPON 862 and the N atoms of
DETDA molecules, similar to the method used by Yarovsky and Evans47 and
Bandyopadhyay et al.20 The modeled crosslinking reaction process is shown in
Figure 2.4. Simultaneous breaking of the CH2-O bonds in the epoxide ends of the
EPON 862 molecules and N-H bonds of the DETDA molecules enable the activated
CH2 ends available to form crosslinks with activated N atoms of the DETDA
molecules. A particular activated N can form a crosslink with the activated CH2 for
any adjacent EPON 862 molecule within a specified cut-off distance. Three
assumptions were made for the crosslinking process:

12
 1) Both primary amines in DETDA were assumed to have the same reactivity

2) The CH2-O and N-H bonds were broken simultaneously (Figure 2.4 A)

3) The newly formed CH2-N and O-H bonds were created simultaneously
(Figure 2.4 B)

Figure 2.4. Crosslinking reaction: (A) Formation of primary C-N bond, leaving a negative charge on
the oxygen and a positive charge on the nitrogen. (B) The negatively charged oxygen abstracts a proton
from the neighboring amine, resulting in an alcohol group and an amine group. The reacted nitrogen is
capable of forming another crosslink by the same reaction, breaking the N-H bond shown in red.

A flow chart describing the modeled crosslink reaction procedure is shown in

Figure 2.5. After selecting an RMS cut-off distance, which affects the ultimate
crosslink density of the MD model and overall energy of the model (larger bond
distances produce higher energy), the RMS distances were calculated between each
pair of N atoms and CH2 united atoms selecting atoms within a specified cut-off
distance. If multiple pairs (for a particular N or CH2) were found within the cut-off,
the closest pair was chosen for the crosslink reaction. Pairs outside the cut-off were
not considered. The covalent bonds corresponding to these crosslink reactions were
inserted into the MD model. The next step of the process was to identify and form
the appropriate secondary bonds between the H atoms from the broken NH2 groups
and the O- atoms of the broken epoxide ends. This bond was formed based on the
closest RMS distances between the O- and H+ atoms from the previous step.

13
Because the newly formed bonds were capable of forming at distances greater than
that of equilibrium, a multistep bond relaxation process was implemented (as shown
in Figure 2.5) to avoid excessively high energies within the system. Newly formed
bonds were run through a minimization and quick NVT simulation with a relaxed
bond force constant for cut-off distances >5 Å.

Figure 2.5. Flowchart describing crosslinking algorithm.

Initial runs of this crosslinking method enabled for high amounts of

crosslinks to be formed for low cut-off distances. The term crosslink density is used
to describe the ratio of the total number of crosslinks that were formed in relation to
the maximum number of crosslinks that could be formed. For example, an epoxy
network having 10 out of 20 crosslinks would have a crosslink density of 50%. For
14
all models described in this work, an initial cut-off of 3.5 Å was incorporated which
achieved an average of ~6% crosslink density in each model. For each iteration
thereafter, the cut-off distance was slightly increased in order to achieve increases in
crosslink density around that original value of ~6% per run. In order to maintain a
realistic crosslinking mechanism, cutoffs did not exceed 6 Å. By not increasing
crosslink density by large amounts per run, the system was able to move more freely
from iteration to iteration. It was discovered early on that by crosslinking large
amounts of polymer (>10% crosslink density per run), the molecules within the
system would become very restricted and successive crosslinking simulations could
not achieve nearly as high of crosslink density. Also, by not having cut-off distances
greater than 6 Å, it significantly reduced the chances of creating a bond through an
aromatic ring. This is known as ‘ring spearing’, and theoretically will not occur in
nature. It is also important to note that as crosslink density increased (>50%
crosslink density), the duration of the NVT simulation for equilibration would
increase. Initially systems were equilibrated for 50ps between runs. For later
crosslink iterations, the NVT equilibration time was increased up to 300ps in order to
better ‘stir’ up the molecules and achieve higher crosslink densities per run. By using
this method, the highest crosslink density achieved was 85% which required 20+
runs of the crosslinking program and required ~5ns of simulation time. Higher
crosslink densities could have been achieved by further increasing the cut-off
distance, though larger cut-off distances become less representative of a ‘realistic’
crosslinking mechanism.

Many different crosslink densities were investigated in this work, as the exact
extent to which a thermoset epoxy crosslinks is unknown. Because of this, it is
important to look at the influence of crosslink density on both molecular behavior
and mechanical properties. Crosslink densities used in each model will be described
in more detail in upcoming chapters.

15
CHAPTER 3: MODELING THE FIBER / MATRIX INTERFACE1

This chapter describes the first models that were created in order to
investigate the interaction between EPON 862-DETDA polymer atoms and a carbon
fiber surface. The effects of crosslink density were taken into consideration. This
initial study utilized a unique method for equilibrating the models to the appropriate
density. All modeling details are highlighted in the following.

3.1 Introduction

The interface between a polymer matrix and carbon reinforcement has long
been a topic of interest. It was recognized early on that a region exists near the
reinforcement surface in which polymer atoms behave differently than bulk.2-4 In
relation to MD, Stevens33 used a course-grained model to simulate interfacial
fracture between a highly crosslinked polymer network and a solid surface.
Mansfield and Theodorou32 investigated the interface between graphite and a glassy
polymer and determined that a 10 Å thick interfacial region existed in the polymer
that was structurally different from that of the bulk polymer. Recently, Chunyu Li et
al.37 used MD simulations to observe the interface of a crosslinked thermoset
polymer in the presence of a graphite surface. Although these studies have provided
valuable information regarding the physical nature of the interfacial region in
composites materials, they have not addressed the influence of crosslink density on
the molecular structure of graphite/epoxy composite interfaces.

The purpose of this study was to determine the effect of crosslink density on
the molecular structure of the fiber/matrix in a graphite/epoxy (EPON 862/DETDA)
composite material. MD models were constructed for a wide range of crosslink
1
Portions of text in this chapter have been reprinted with permission from C.M. Hadden, B.D. Jensen,
A. Bandyopadhyay, G.M.Odegard, A.Koo, R.Liang Composite Science and Technology 2013, 76, 92-
99. © Elsevier Ltd. Please refer Appendix B for copyright information.
16
densities and the molecular structure of the interface corresponding to each crosslink
density was determined. The simulations predicted a polymer molecular structure at
the interface that is different than that of the bulk polymer and is dependent on the
crosslink density.

3.2 Modeling Procedures

The initial un-crosslinked polymer model began as the 2:1 stoichiometric

ratio of the EPON 862 monomer and the DETDA hardener outlined in section 2.2
(page 9). The initial 2:1 structure was then implemented into a program that
randomly rotated the molecular coordinates, and then translated these randomized
coordinates along the x, y, and z axes to fill a specified volume. The resulting
structure consisted of many randomly oriented clusters of the small 2:1 ratio stacked
loosely together in a manner much like that of a simple cubic crystal structure, as
seen in figure 3.1. The 2:1 coordinates file was multiplied into a low density
structure consisting of 15,936 polymer atoms with a ratio of EPON 862:DETDA of
394:192 and a simulation box size of 64 x 66 x 115 Å, equating to a density of ~0.5
g/cc. The larger polymer model was then allowed to equilibrate into a uniform low-
density liquid using a slow relaxation procedure performed over a cycle of 6 MM
and 3 MD simulations. All MD simulations were conducted with the NVT ensemble
for 100ps at 300K. This equilibration was done before the addition of a simulated
graphite surface.

17
Figure 3.1 Resulting structure after implementing the ‘randomly rotate and translate’ program for
multiplying the initial 2:1 ratio structure.

For the creation of the graphite surface, a program was written to create
sheets of carbon based on the aromatic pattern for which graphite surfaces are well
known. By assigning atom IDs in a specific pattern, sections of code were easily
written to identify the appropriate bonds, angles, and dihedrals for a given sheet of
aromatic carbon geometry. The simulated graphite surface was constructed from 3
sheets of stacked graphene, each sheet containing 1,728 carbon atoms for a total of
5,184 atoms. The graphene sheets were oriented along the x-y plane, with periodic
boundary conditions in the x and y-direction, and had an interlayer spacing of 3.35Å.
Sheets were initially planar and contained no terminal hydrogens, surface
imperfections, nor surface treatments. The graphite structures were relaxed using a
series of MM and MD simulations, similar to that described above for the polymer
structure. This initial equilibration step was performed without the presence of the
polymer molecules. While equilibrating the graphite, the z-direction box coordinate
was chosen to implement interlayer spacing for periodic boundary conditions. Thus
the top surface was influenced by the bottom surface and visa-versa, representing
18
many layers of bulk graphite. Initial and relaxed graphite structures are shown in
Figure 3.2. Upon relaxation, the graphite structure was positioned to line up with the
64 x 66 Å dimension of the polymer structure.

Figure 3.2 a) Initial planar graphite structure and b) Graphite structure after equilibration

For the next step in the procedure it was desired to simulate the non-
crosslinked, low-density epoxy liquid resin in a parallelepiped simulation box
surrounded on the bottom side (lower z-coordinate) by the graphite structure, on the
opposite side (upper z-coordinate) by the bulk epoxy resin, and on the remaining
four sides (x and y) by replicate images of the low-density liquid resin. This is shown
on the left side of Figure 3.3 where the simulated polymer (labeled as “polymer”) is
situated above the graphite molecules along the z-axis and below a bulk polymer
molecular structure. While the graphite and low-density polymer structures were
established as described above, the bulk polymer structure is a fully crosslinked
model of 76% crosslinked EPON 862/DETDA resin with a bulk density of 1.2 g/cc
and 17,928 united atoms that was obtained previously by fellow researcher, A.
Bandyophadyay.20 Periodic boundary conditions were assumed along the x- and y-
axes. Therefore, the molecular model effectively simulated an infinitely large, flat
interface. Although fiber surfaces are round, the radius of curvature at the molecular
level is large enough to effectively model it as a flat surface. An initial gap of 3.35 Å
(interlayer spacing of graphite) was placed between the graphite and polymer. The
entire molecular model was then allowed to partially equilibrate using two MM
simulations and one 100 ps MD simulation, so that molecules from each individual
19
structure (polymer, graphite, bulk) could react to one another The entire model
contained a total of 39,048 united atoms.

Following this initial equilibration step, the atoms of the bulk polymer and
the graphite sheets were fixed (no movement) and only non-bonded interactions
between them and the polymer were considered for the remaining simulation
processes. This step was performed to reduce computation time and to focus the
simulations on the influence of the presence of the graphite sheet on the adjacent
polymer molecular structure.

During the next step of establishing the composite interface molecular model,
the simulated polymer molecules were condensed to achieve the desired polymer
density. The non-bonded van der Waals parameters were initially scaled down to
prevent large increases in the energy of the system. In order to reach the appropriate
polymer density, the atoms of the bulk polymer were displaced in small incremental
amounts using the following equation,

zinew zicurrent  S  zicurrent  z bottom  (3.1)

where zinew is the new z-coordinate for polymer atom i, zicurrent is the current z-
coordinate of polymer atom i, zbottom is the lower-most polymer atom z-coordinate,
and S is a scaling factor for the amount of displacement desired. For this study, S
varied between 0.02 and 0.05, becoming smaller as the system became closer to the
target density and the non-bonded van der Waals parameters were gradually scaled
up to normal values. Under this constraint, the polymer atom with the lowest z-
coordinate value (nearest to the graphite) was not displaced and the polymer atom
with the highest z-coordinate (nearest to the bulk polymer) was displaced the most
for each incremental step. The bulk polymer molecules were uniformly displaced
the exact amount as the polymer atom with the highest z-coordinate, thereby
continually condensing the polymer atoms. Each displacement was followed by a
MM minimization and a 50ps MD simulation, followed lastly by another MM
minimization. If a particular displacement step resulted in total energy values that

20
were relatively high, the previous equilibration step was repeated before further
condensation steps were taken. This process continued until achieving a polymer
density of 1.16g/cc, totaling over 1 ns of MD simulation time. The resulting
molecular model is shown on the right side of Figure 3.3. It should be noted here that
this initial condensation procedure was implanted before the discovery of the
fix/deform command in the LAMMPS software, which allows the user to deform the
simulation box along any specified axis during an MD simulation. The fix/deform
command was implemented in later models.

Figure 3.3 Combining the three separate MD models (left) to form the complete MD model of the
composite interface (right).

After condensing the model to the desired density, crosslinking was

performed using the method described in section 2.3 (page 12). A total of six
molecular systems were established, each having a unique crosslink density. The
chosen crosslink densities were 57%, 65%, 70%, 75%, 80%, and 85%. After
crosslinking to the desired density, each structure was allowed to equilibrate using 3
MM minimizations and 2 MD simulations of 1 nanosecond. The density of formed
21
crosslink atoms (C-N and O-H) as a function of the z-axis can be seen in Figure 3.4,
where the origin of the z-axis lies in the center of the top graphene sheet, as shown
on the right side of Figure 3.3. The figure shows good dispersion of crosslinks
throughout the polymer structure with unique profiles for each crosslink density.
Therefore, the crosslink distributions are relatively uniform and independent of
crosslink density. All density calculations were performed using the fix/ave spatial
command in the LAMMPS software.12

Figure 3.4 Crosslink atom mass density along the z-axis for each crosslinked structure. Z-axis origin
is taken from the center of mass of the closest graphene sheet.

3.3 Results and Discussion

After the models were equilibrated, they were used to examine the effects of
crosslink density on the overall polymer mass density at the polymer-graphite

22
interface, the internal stress concentrations in the polymer at the interface, and the
different potential energy components in the polymer.

Polymer Mass Density Characteristics

Mass densities were computed using the “fix/ave” spatial command in

LAMMPS12, which sums the per-atom densities and averages them for slices of a
specified thickness along a specified axis. A slice thickness of 0.2 Å stacking along
the z-axis was chosen to observe the structural change moving away from graphite
surface oriented in the x-y plane. Figure 3.5 shows the overall mass density profile for
the three graphite sheets and polymer molecules at the interface. The finite width of
the graphite sheet peaks is due to their slightly non-planar shape (Figure 2). The mass
density of the polymer near the surface differs from that in the bulk with a profile
similar to those reported for polymer systems with graphite27,33 , carbon nanotube25
and nanoparticle29,30 reinforcement. Common interfacial characteristics include
(moving along the positive z-axis) an initial peak above bulk mass density, followed
by a trough below bulk density, ending with a small peak above bulk density before
leveling off.

Figure 3.5 Mass density profile along Z-axis of simulation box for graphite and polymer molecules.
Z-axis zer value is taken as the center of mass for the outermost graphite layer. Polymer bulk density
is shown by the red line (~1.16g/cc).
23
Figure 3.6 Interfacial mass density profile along Z-axis of simulation box (zoomed in from Figure
3.5)

Figure 3.6 shows the same data shown in Figure 3.5, but focused on the
interfacial region within 15 Å of the graphene sheet. The fluctuation in mass density
is observed for about 10 Å from the center of the nearest graphite sheet. Therefore, the
effective surface thickness of a polymer in the vicinity of graphite is about 10 Å. Figure
3.6 also demonstrates that the surface effects are reduced with increasing crosslink
densities. As the polymer approaches higher crosslink densities, the mass density peak
amplitudes reduce and the structure more closely resembles the bulk mass density
distribution. This behavior is most likely due to the tendency of the crosslinks to hold
the network together in a more spatially consistent manner.

The predictions shown in Figures 3.5 and 3.6 are supported by the transmission
electron microscopy (TEM) image shown in Figure 3.7. The image shows a portion of
an open-ended carbon nanofiber reinforcing a crosslinked epoxy matrix. The surface
of the carbon nanofiber and the structure of the epoxy system are similar to those
modeled herein. The bulk epoxy matrix is amorphous; however, a small band of
24
structured epoxy can be seen on the surface of the nanofiber, as indicated in the figure.
Although the thickness of this interface appears to vary in the image, the magnitude of
the thickness is close to the predicted value, 10 Å. This favorable comparison serves
as a validation of the modeling strategy described herein.

Figure 3.7 TEM image of carbon nanofiber embedded in a crosslinked epoxy resin. Image courtesy of
Z. Liang et. al, High-Performance Materials Institute, Florida State University

25
 Figure 3.8 shows the mass density profiles for the EPON 862 and DETDA
molecules for the non-crosslinked and 80% crosslinked systems. In general, the
densities of EPON 862 are larger than those of DETDA, which is mostly because there
are two EPON 862 molecules for every DETDA molecule. For both levels of
crosslinking, larger concentrations of EPON are observed near the surface. Large
concentrations of DETDA are only present at the surface for the non-crosslinked
system. Therefore, before crosslinking occurs, both molecules are concentrated at the
surface. During the process of crosslinking, the concentrated DETDA molecules are
pulled away from the surface so that the DETDA mass density profile is relatively
uniform.

Figure 3.8 EPON 862 and DETDA molecular mass densities along simulation box Z-axis for Non-
crosslinked and 80% crosslinked structures.

Polymer Internal Stress Characteristics

Stress components were computed using the LAMMPS “fix/ave” spatial

command in a similar manner as described above for the mass density profiles. Each

26
normal and shear stress component was computed for each individual atom and
averaged in slices of 0.1 Å along the z-axis of the simulation box. The individual
stress components were then used to calculate the Von Mises stress using the following
equation:

ଵ
ߪ௏ெ =
ξଶ
ඥ(ߪ௫ െ ߪ௬ )ଶ + (ߪ௬ െ ߪ௭ )ଶ + (ߪ௭ െ ߪ௫ )ଶ + 6(߬௫௬
ଶ + ߬ ଶ + ߬ ଶ ) (3.2)
௬௭ ௭௫

The Von Mises stress along the z-axis for a range of crosslink densities is
shown in Figure 3.9. The stress peak nearest to the origin is due to stresses in the
graphite sheet. Residual stresses are observed for each of the crosslink densities with
the largest values at the interface. Despite the large amount of scatter in the data, there
appears to be no influence of the crosslink density on the location and magnitude of
residual stresses. The largest amount of stress in the polymer is considerably low for
all crosslink densities, well below the expected tensile strength of EPON-862, which
is 70-95 MPa at room temperature.48 The presence of a small concentration of residual
stress near the interface indicates that upon significant loading of an epoxy/graphite
composite, failure of the resin matrix would likely occur first near the fiber/matrix
interface.

0.5
Non-crosslinked
57% crosslinked
0.4 65% crosslinked
70% crosslinked
Von Mises Stress (MPa)

75% crosslinked
0.3 80% crosslinked
85% crosslinked

0.2

0.1

0
0 10 20 30 40 50
z-zxis (Angstroms)

Figure 3.9 Spatially averaged Von Mises stress for polymer and graphite atoms along Z-axis.

27
Polymer Potential Energy Characterisitics

Potential energy values for bonds, angles, dihedrals, and non-bonded

interactions were computed with LAMMPS using the same fix/ave spatial command
previously described. The bin size was 0.1 Å along the z-axis. Figure 3.10 shows the
total potential energy (the sum of the individual energy terms) for each system. There
appears to be an increased level of potential energy near the fiber/matrix interface,
which is consistent with the increase in stress at the interface shown in Figure 11. Also,
the potential energy generally decreases with increasing crosslink density. The largest
contributing factor to the total potential energy was found to be angle energy, which
decreases with increasing crosslink density. This is likely due to the elimination of the
high-energy epoxide rings on the EPON 862 molecule during the crosslinking
reaction. The decreases in potential energy with increasing crosslink density shown
in Figure 12 indicate that crosslinking is energetically favorable for the modeled
system and force field.

Figure 3.10 Spatially averaged per-atom potential energy. Potential energy values are the sum of bond,
angle, dihedral, and non-bonded van der Waal’s energies.

28
3.4 Conclusions

In this chapter, the influence of crosslink density on the molecular structure of

the graphite fiber/epoxy matrix interface was examined using MD techniques. The
modeled epoxy resin was an EPON 862/DETDA system. It was determined that the
mass density of the polymer within 10 Å of the center of the nearest graphene sheet
was perturbed from the bulk level of 1.17 g/cc. This observation agrees with
predictions in the literature and our own TEM images (Figure 3.8). Although the
amount of crosslink density did not influence this effective surface thickness of the
epoxy near the graphite interface, increasing levels of crosslinking reduced the
magnitude of mass density fluctuations within the effective surface thickness. Analysis
of the mass density of polymer molecules also reveals that before crosslinking, the
presence of surface effects is due to both the perturbed mass densities of the EPON
862 monomer and DETDA hardener molecules in the surface region. Local Von Mises
stresses were determined at the interface of the composite, and the MD simulations
demonstrated that the internal stresses were slightly higher in the surface region than
in the bulk epoxy, albeit at levels far below the ultimate strength of the neat resin. The
internal stresses near the interface are not influenced by crosslink density. Finally, the
MD simulations predict elevated levels of molecular potential energy in the surface
region of the polymer, with the potential energy magnitude decreasing with increasing
crosslink levels. It can be said, with confidence, that the methods described in this
chapter are sufficient for producing accurate models of EPON 862-DETDA epoxy in
the presence of a graphite surface.

29
 CHAPTER 4: MULTI-SCALE MODELING OF GRAPHENE
DOPED EPOXY 1

Chapter 4 utlizes MD to investigate the effects of graphene thickness along

with varying cross-link density on the molecular scale mechanical properties. The
properties derived from MD are then implemented into a Micromechanical Analysis
Code with the Generalized Method of Cells (MAC/GMC) software for comparison
to larger scale experimental results for materials similar to the models.

4.1 Introduction

The use of graphene nano-platelets (GNPs) within the epoxy matrix has been
shown to enhance both mechanical and electrical properties within the composite.
The effects of GNPs on the composite is governed by the amount of GNPs added to
the matrix component and the dispersion of the GNPs within the matrix9-10,49-50. The
interface between GNP and surrounding epoxy matrix is difficult to characterize with
experimental techniques, and MD can provide valuable insight to the material
behavior at the molecular level. MD can also provide information for GNPs of
different sizes (atomic thicknesses), as GNPs vary significantly in size depending on
fabrication techniques.

Recently, Chunyu Li. et. al. utilized MD techniques to model a multi-layer

graphene epoxy composites. Graphene orientation was investigated, though the
number of atomic layers and crosslink density of the polymer remained constant.
Mechanical properties were also derived, though no experimental comparison was
incorporated37. Jing-Suei Gao et. al. investigate the mechanical properties of poly-
imide polymer near a graphene sheet of 5 atomic layers thick. Interfacial effects were

1
The information from this chapter will be considered for publication in the near future.
30
investigated and mechanical properties were derived. Again, the results did not allow
for comparison to larger scale experimental data.36 Although these studies have
given valuable information regarding the physical nature of the interfacial region,
including an estimate of the mechanical properties on a molecular level, there has
been little effort to implement this molecular data into a larger scale model that will
enable comparison to experimental results for similar systems.

The purpose of this study was to use mechanical data provided by MD for
epoxy/graphene systems, and implement this data into a larger scale model using the
principals derived from the Generalized Method of Cells (GMC).51 MD models were
constructed for a wide range of epoxy crosslink densities (curing ranges) with
graphene structures of varying atomic thickness from a single atomic layer up to 4
layers thick. Mechanical properties were extracted from each model and then
inserted into Micromechanics Analysis Code with Generarlized Method of Cells
(MAC/GMC) software51-53 to calculate properties for an epoxy/GNP composite with
varying GNP concentration and dispersion at a fixed polymer cross-link density
(explained in detail in section 4.3). The results for a GNP doped epoxy matrix for a
specified GNP concentration and dispersion can be re-introduced to the MAC/GMC
software in the presence of a carbon fiber to yield properties for a GNP doped
epoxy/carbon fiber composite. MAC/GMC results for the GNP doped epoxy matrix
calculation, along with that for a GNP doped epoxy/carbon fiber composite are in
good agreement with experimental data. Composite properties increase with an
increase in GNP concentration. Results are best for single layered GNPs, which also
represent the highest amount of dispersion in the GMC models.

4.2 Modeling Procedures

The initial un-crosslinked polymer model began as the 2:1 stoichiometric

ratio of the EPON 862 monomer and the DETDA hardener outlined in section 2.2
(page 9). The 2:1 molecular model was subjected to four MM minimizations and
31
three 100 ps MD simulations. MM simulations utilized the conjugate gradient
stopping criterion, and MD simulations were performed using the NVT ensemble at
300K.12 This process minimized internal forces and thus reduced internal residual
stresses which were created from the initial construction of bonds, bond angles, and
bond dihedrals. After the structure stabilized to a relatively low energy value, the
initial 2:1 stoichiometric structure was replicated, and the replicated models were
randomly rotated and then translated along the x, y, and z axes and combined into a
much larger structure with an EPON 862:DETDA ratio of 250:125, containing
10,375 total united atoms. Therefore, the resulting system consisted of 250
randomly oriented clusters of the small 2:1 ratio cluster stacked loosely together in a
manner much like that of a simple cubic crystal structure. This larger polymer model
was then mirrored about a graphene structure positioned in the x-y plane central to
the z-axis. As a result, each system contains a 500:250 ratio of EPON 862:DETDA
totaling 20,750 polymer atoms. The centralized graphene structures varied in atomic
thickness from 1 atomic layer up to 4 layers thick, each layer containing 4200 carbon
atoms. The largest system, comprised of a 4 layer graphene sheet, contained 37,550
total atoms and the initial box size was 101 x 104 x 210 Å. All models employed 3D
periodic boundary conditions, which is necessary for deformation simulations. The
initial box size produced a polymer density approximately equal to half of a fully
cured solid EPON 862 epoxy (~0.5 g/cc in all 4 systems). Unlike the previous model,
the mirrored epoxy structures and the graphene structures were all combined in one
step and allowed to equilibrate in the presence of one another for 200ps. In order to
achieve the desired polymer density of 1.17 g/cc, the 4 separate models were
subjected to 12 cycles of deformation along the z-axis. Each cycle included a MM
followed by a 100ps MD NVT simulation in which the z coordinate was reduced in
equal amounts from both the positive and negative z-coordinate boundaries using the
LAMMPS fix/deform tool. A Nose/Hoover thermostat and barostat was
implemented for temperature and pressure control, respectively.46 The amount of
deformation decreased with each cycle as the models became closer to the desired
density. This was done to avoid large energy increases to the system by packing the

32
molecules together too quickly. This entire densification process was done over a
total of 1.2 ns for each of the 4 systems. The final z-coordinate boundary enabled for
polymer atoms to extend ~13 Å from the graphene surface, to ensure that the
interfacial region was fully captured, and show a minimal influence from the bulk
polymer characteristics during deformation. A visual representation for the
densification process can been seen for the single graphene sheet system in Figure
4.1. Fully equilibrated, non crosslinked, structures for all 4 systems can be seen in
Figure 4.2.

Figure 4.1 Molecular structures for single layer graphene sheet model highlighting the densification
process.

33
Figure 4.2 Equilibrated models for varying number of graphene layers before crosslinking. Polymer
density ~1.17 g/cc for each model.

After condensing the models to the desired density, crosslinking was

performed as previously described in section 2.3. For this study, multiple crosslink
densities were considered in conjunction with varying graphene thickness. The chosen
crosslink densities were 65%, 70%, 75%, and 80%. Figure 4.3 shows the density
profile along the z-axis for the 80% crosslinked models. Again, there is good
dispersion of crosslinks throughout the models with unique profiles for each model
with a varying number of graphene atomic layers. Crosslink density profiles for 65%,
70%, and 75% produce similar results.

34
Figure 4.3 Crosslink dispersion as a function of z-axis coordinate for 80% crosslinked systems.

Mechanical Deformation

Molecular models were assumed to be transversely isotropic in the x and y

directions. In order to develop Stiffness matrices for each of the 16 systems, models
were deformed in both tension and shear. Tensile deformations of 5% were performed
separately along the x, y, and z axes over the period of 1ns for each model. Tensile
deformations aloud for Poisson’s contractions in the transverse directions for
calculation of Poisson’s ratio. In order to accomplish this, it is necessary to utilize the
NPT ensemble to allow for volume adjustments. Shear deformations of 5% were
performed separately along the x-y, y-z, and x-z plane over the period of 1ns for each
model. Shear deformations do not allow for contraction, and thus incorporate the NVT
ensemble. A visual representation of tensile deformation along the z-axis, and shear
deformation in the x-z direction is shown in Figure 4.4. Resulting stress vs. strain
curves for tensile deformation along the x-axis and shear in the xy-plane are shown
for the 80% crosslinked models in Figure 4.5. Data shows very little scatter, indicating
an acceptable strain-rate. All mechanical data derived from MD models can be seen in
Table 4.2 shown later in this chapter.

35
Figure 4.4 Simulated deformation characteristics for tension along the z-axis (left) and shear in the x-
z direction (right)

Figure 4.5 Tensile and Shear stress vs strain curves for varying number of graphene layers in 80%
crosslinked structures.

36
Micromechanics

Larger scale models were developed using MAC/GMC software created by

Brett Bednarcyk et. al.51-53 The software utilizes the concepts of the Generalized
Method of Cells (GMC) from micromechanics.51 GMC estimates macro-scale
material properties by investigating a very small representative unit cell (RUC) that
can be multiplied in a periodic array into a larger scale composite. The RUC
commonly contains volumes representing a polymer matrix and a carbon
reinforcement, where the material properties from each component determine that of
the entire composite.

GNP Doped Epoxy Model

In order to develop a model for a GNP doped EPON 862 epoxy, the mechanical
properties derived from each MD model were implemented into a small volume
portion of the RUC. The surrounding volume in the RUC represented that of fully
cured bulk EPON 862-DETDA epoxy. This 3D RUC was then randomly rotated and
translated in the x, y, and z directions to create a large cube with many randomly
oriented graphene structures, as seen in Figure 4.6. The GMC model allows for control
of the graphene volume concentration for the overall composite, and the randomly
oriented RUC was multiplied extensively to ensure a good representation of an
isotropic GNP doped epoxy matrix. As the number of graphene atomic layers increases
in the MD models, so does the graphene volume fraction for each RUC. Therefore, if
the graphene volume within the RUC is kept constant, the model containing data from
1 atomic layer of graphene would have to be nearly 3 times as large as that for the
model containing data from the 4 atomic layer graphene structure in order to achieve
the same graphene concentration in the GMC models. Also, the graphene
concentration in the GMC model is governed by the graphene volume fraction in the
MD model. For the MD model with a single atomic layer of graphene the volume
fraction is 0.11, the GMC model will not be able simulate a higher composite graphene
concentration than this. When thinking of the amount of GNPs within a composite, a
single atomic layer of graphene can be thought to represent excellent dispersion of

37
GNPs within the epoxy matrix, while the 4 atomic layer graphene structures will
represent agglomerated GNPs within the matrix component. Volume fractions of
graphene can be seen for each of the MD structures in Table 4.2.

Figure 4.6 Schematic representation for implementing MD results into MACGMC software to obtain
data for GNP doped epoxy and GNP epoxy/Carbon Fiber composites

GNP Epoxy / Carbon Fiber Composite

After developing micromechanics models for isotropic GNP doped epoxy

specimens, the data was re-implemented into MAC/GMC in the presence of a carbon
fiber. This step is shown on the right side of Figure 4.6. MAC/GMC has the ability to
consider many different unit cell architectures for a carbon fiber. For this study, the
fiber architecture ID (ARCHID) was chosen as a 26x26 circular array seen Figure 4.7.
This figure represents the RUC used for composite data calculations. The outer portion
(green) represents the GNP doped epoxy, and utilized input parameters derived from
the techniques described in the previous section. The circle in the center of the RUC
is the chosen geometry to represent the carbon fiber. Input parameters for the carbon
fiber were chosen to accurately represent the fibers used in related experiments, and
are shown in Table 4.110,54

38
 Table 4.1 MAC/GMC Input Parameters for HexTow® AS4-GP/3K Carbon Fibers

Figure 4.7 MAC/GMC RUC for GNP doped epoxy / carbon fiber model. (ARCHID = 13 for
MAC/GMC version 4.0)

4.2 Results and Discussion

MD models were used to examine material behavior at the interface. Upon

validating each model with past simulation results, MD was used to deform the
simulation box and extract specific material properties. These properties were then
implemented in to MAC/GMC software to estimate macroscopic material properties
for both a graphene doped epoxy matrix with varying graphene concentrations, and a
graphene doped epoxy matrix with unidirectional carbon fibers with volume
fractions similar to experimental data.

39
Molecular Dynamics

Polymer mass density profiles follow a similar trend to previous works,

showing peak values near the graphene surface then fluctuating and eventually
becoming stable at ~10 Å from the surface. Figure 4.8 shows the density vs. z-axis
coordinate for varying number of graphene layers and crosslink densities. The larger
spikes in the center of the density profile are from the graphene carbon atoms. These
peaks are not perfectly symmetrical because the graphene is able to wave and contort
during MD. Exact values for overall polymer density at the interface are given in Table
4.2. These values are ~.04 g/cc less than the target value of 1.17 g/cc, likely caused by
the final 1.2ns NPT relaxation simulation described above.

Figure 4.8 Atom density profile along z-axis for a)1 layer b) 2 layers c) 3 layers and d) 4 layers of
graphene

Mechanical deformation results for each system can be seen in Table 4.2.
Young’s modulus increases with increasing layers of graphene but shows little
influence from cross-link density. Deformation simulations in the x-y plane produced
stress vs. strain curves with very little scatter due to the fact that most of the force
40
was placed upon the graphene structures. A visual representation of this can be seen
in Figure 4.9, showing atomic stress values for tensile deformation parallel to the
graphene surface. Shear values for x-z and y-z plane deformation along with
Poisson’s ratios for Z-axis deformation show a lot of scatter because van der waals
attraction between polymer and graphene atoms allows them to easily slide past one
another. Therefore, shear values for Gxz = Gyz were very close to zero and were
given a value of 1 MPa for all models. In order to say that the models were truly
transversely isotropic (i.e. Exx = Eyy, Vxy = Vyx, etc.), which requires less input for
the MAC/GMC software, all applicable isotropic values were taken to be the average
between the two.

Table 4.2 Resulting data from MD simulations

Figure 4.9 Atomic stress before and after tensile deformation along the y-axis.

41
Micromechanics using MAC/GMC software

The mechanical properties derived from MD were implanted into the

MACGMC software as a representation of a GNP with varying layers surrounded by
interfacial polymer atoms. The amount of GNP that would contribute to the
micromechanics calculations was controlled by volume fraction of graphene. Again,
this calculation would incorporate the graphene volume fractions from MD for varying
number of layers. The elastic modulus for bulk EPON 862 at room temperature was
taken to be 2.72 GPa10. Using the micromechanics software, representations of GNP
doped EPON 862 were created with up to 5% volume fraction of GNPs for varying
number of graphene layers and crosslink densities. Results are shown in Figure 4.10,
and indicate that GNPs with fewer layers will have a greater effect on the mechanical
properties of this particular system. The single layer graphene systems were also the
only systems to show any influence from increasing number of crosslinks (curing),
likely because the MACGMC calculations for GNP volume fraction would need to
incorporate more contributions from MD models for a single layer of graphene when
compared to MD contributions for multi-layer models. Because crosslink density had
little effect on the macroscopic properties predicted by MAC/GMC, only the data for
80% crosslinked structures would be considered for the construction of the GNP doped
epoxy/carbon fiber micromechanics models described later.

42
Figure 4.10 MAC/GMC data for addition of GNPs with carrying number of layers and interfacial
polymer crosslink density.

The experimental data used for comparison was taken from Dr. J. King et al.10
where tensile test specimens were constructed with EPON 862 / DETDA epoxy and
xGnP®-C-300 GNPs from XG sciences.54 The average GNP thickness was 2nm before
mixing, which would equate to 7-9 layers of graphene. The mixing process may also
further break down the GNPs and reduce the overall particle thickness. Experimental
results agree well with the data generated from the MD models with 3 or more layers
of graphene, as seen in Figure 4.10. Results for the 1 and 2 layer GNPs produce
significant increases in Elastic Modulus, though fabrication techniques for producing
GNPs with a thicknesses less than 1nm have not yet been developed.

Results for the MAC/GMC models considering GNP doped epoxy/carbon fiber
composites are shown in Figure 4.11, and show a slight increase in modulus values for
increasing GNP volume fractions. As expected, results for GNPs of a single atomic
layer produce the greatest increases in modulus values. Experimental results for
comparison were taken from Dr. King et al. for test specimens constructed of
HexTow® AS4-GP/3K carbon fibers with the same xGnP®-C-300 GNPs used in the
43
epoxy/GNP specimens previously mentioned. Experimental results are shown in Table
4.3, and indicate a decrease in tensile modulus with increasing GNP volume fractions.
Though the model data and the experimental data observe different trends, the model
data falls within the standard deviation of the experimental results for all GNP volume
fractions, as shown in Figure 4.12. Transverse modulus results from MAC/GMC
observe a similar trend as tensile modulus, increasing with higher GNP volume
fractions and the largest increases are seen in the single layer model. Transverse
modulus values are shown in Figure 4.13.

Figure 4.11 MAC/GMC tensile modulus for GNP doped epoxy/carbon fiber models for varying GNP
volume fractions and atomic thicknesses.

44
Table 4.3 Experimental results for epoxy/carbon fiber composites with varying GNP volume fractions

Figure 4.12 Experimental and MAC/GMC tensile modulus for GNP doped epoxy/carbon fiber models
for varying GNP volume fractions and atomic thicknesses.

45
Figure 4.13 MAC/GMC transverse modulus for GNP doped epoxy/carbon fiber models for varying
GNP volume fractions and atomic thicknesses.

4.4 Conclusions

In this study, MD techniques were utilized to determine the influence of crosslink

density and GNP atomic thickness on molecular level mechanical properties. The
models were designed to focus on the interfacial polymer surrounding an embedded
graphene structure with different numbers of atomic layers. Models show an
interfacial region of ~10 Å in which the mass density profile behaves differently than
the bulk polymer, further validating previous works. Mechanical properties derived
from MD models such as tensile modulus and shear modulus increase with
increasing number of graphene atomic layers. Varying crosslink density shows no
significant effect on mechanical properties for the range investigated by MD models.
Micromechanics models were created for GNP doped epoxy specimens using
MAC/GMC software, with inputs from MD results. The micromechanics analysis
was designed to create a fully isotropic representation of a GNP doped epoxy

46
specimen with GNPs of varying atomic thickness and GNP volume fraction. Results
show increases in Elastic Modulus for increasing GNP volume fractions, with the
greatest increases in modulus coming from GNPs of a single atomic layer and
decreasing as more layers are considered. Results from MAC/GMC for GNP doped
epoxy specimens are in good agreement with experimental data, especially for
modelled GNPs of 3 or more atomic layers. MAC/GMC models were next created
for GNP doped epoxy/carbon fiber composites. These composite models show slight
increases in both tensile and transverse modulus values with increasing GNP volume
fractions in the surrounding matrix epoxy. Again, GNPs of a single atomic layer
produced the most significant increase. Experimental comparison shows the
composite models to fall within the Standard Deviation for all sampled GNP volume
fractions and further validate this method for utilizing MD data to develop larger
scale models of composite materials.

CHAPTER 5: MULTI-SCALE MODELING OF GRAPHENE

OXIDE / EPOXY COMPOSITES2

The fabrication of carbon fibers and graphene nanoplatelets (GNPs) will

often leave oxygen species bonded to the surface of the carbon structures. In this
chapter, different oxygen species are attached in different amounts to a single layer
graphene surface and a 4 layer graphene surface. The results are compared to
previous data for pristine graphite structures.

5.1 Introduction

The surface of graphite structures often contains small amounts of oxygen

and even smaller amounts of nitrogen that are deposited during the fabrication
process.55 These groups can influence the behavior of the interfacial region between

2
The information from this chapter will be considered for publication in the near future.
47
the surface of the graphite structure and the surrounding epoxy matrix. There are
many accepted models for the structure of graphene oxide, as shown in Figure 5.1.56
Experimental techniques such as X-ray photoelectron spectroscopy (XPS) can help
reveal the exact oxygen species present, along with a rough estimate of their atomic
composition for the tested graphite structure.57 The use of molecular dynamics (MD),
coupled with knowledge gained from experimental data, can be a powerful tool in
discovering the effects oxygen on the surface of a carbon structure. MD is able to
provide insight into: the effects of adding specific oxygen species, the effects of
varying the amount of these species on the graphite surface, the likely pattern or
confirmation of which these oxygen species will attach to the surface of a graphite
structure. All of which are still topics of relative uncertainty and or unpredictability.

Figure 5.1 Widely accepted models for the structure of graphene oxide.

48
 The use of MD in relation to graphite oxide has been explored by several
researchers.58-64 Early models by Allington et. al. incorporate regularly spaced
oxygen species (OH, CO, COOH) on a graphitic surface and report their interaction
energies with small molecules58,59 I.S. Awan et. al. utilized MD to construct models
of oxidized carbon fiber surfaces, varying concentrations of various oxygen species
and reporting the changes is surface energy. These studies have provided valuable
information regarding the molecular structure of oxidized graphite, though there has
been little effort to investigate the effects on mechanical properties.

The purpose of this study was to create models of oxidized graphene surfaces
surrounded on both sides by EPON 862/DETDA epoxy. The data provided from
referenced MD research would become the basis for deciding which oxygen species
to add, the concentration of oxygen species upon the graphite surface, and the
bonding arrangement on the graphite surfaces. Models were created for graphene
oxide with the addition of epoxide oxygen species, hydroxyl oxygen species, and a
combination of epoxide and hydroxyl oxygen species. These three models were
created for graphene structures of 1 and 4 atomic layers. The amount of oxygen
present on the surface was also varied for amounts of 12% and 24% oxygen. The
influences of these 3 variables are investigated within the MD models. Mechanical
properties are then derived from the MD models and compared to results for pristine
graphene models from Chapter 4. Results show oxygen groups to cause a decline in
modulus values in all cases. The addition of epoxide oxygen species show the least
amount of reduction when compared to pristine graphene structures.

5.2 Modeling Procedures

The molecular models described in this chapter were created nearly exactly
the same as those from Chapter 4. The initial 2:1 Stoichiometric ratio of EPON
862:DETDA was replicated into a 250:125 ratio, and this structure was mirrored

49
about a central graphene oxide structure. Therefore, the final ration of
EPON:DETDA was 500:250 totaling 20,750 polymer atoms. Each layer of graphene
was again composed of 4200 carbon atoms, and the box size was kept the same. The
focus of this chapter was the effect of oxidizing the graphene surface. In order to
save time, graphene structures were created composed of only a single atomic layer
and of 4 atomic layers. This would enable comparison to the optimal result from the
previous chapter (single layer graphene) and the most realistic result from the
previous chapter (4 layer graphene). The oxygen species chosen for this project were
based on a study performed by D. W. Boukhvalov et. al., in which Density
Functional Theory (DFT) was used in conjunction with the widely accepted models
of graphite oxide (Figure 5.1) to identify the most stable attachment of oxygen
species to graphite structures.61 From this research, The three most stable
configurations for the attachment of oxygen species to a graphite surface were
chosen to be implemented into the MD model. They can be seen in Figure 5.2.
Models were created with epoxide molecules only (Fig. 5.2a), hydroxyl molecules
only (Fig 5.2b), and a combination of epoxide and hydroxyl molecules (Figure 5.2c).
The lowest energy configuration, as found by DFT, was the epoxide only model. The
next lowest energy was that for the hydroxyl only model and the highest energy
confirmation was the combination model.61 Each of these oxygen species would be
attached to the single sheet models, along with the outer most layer of the 4 sheet
models, in a regular pattern in varying amounts. A total of 12 models were created
varying in: number of graphene layers (1 and 4 layered structures), the attached
oxygen species (epoxide, hydroxyl, and a combination), and the amount of each
oxygen species. The amounts of oxygen are described in this research as a
percentage, and correspond to the number of carbon atoms that they are attached to.
If oxygen species were attached to 420 of the 4200 carbon atoms, that would equate
to 10%. The amounts chosen for this work were 12% and 24%, which are in good
agreement with referenced data.55-64 Figure 5.3 shows the fully equilibrated graphene
oxide structures of both the single and 4 atomic layers for the 3 possible
combinations covering 12% of the carbon atoms. As shown in Figure 5.2, oxygen

50
species protrude from both sides of the graphene surface. This is difficult to see from
Figure 5.3, but the focus of the figure was the regular pattern in which oxygen
species were bonded.

Figure 5.2 The most stable confirmations of graphene oxide as found by DFT for a) epoxide oxygen
only, b) hydroxyl molecules only, and c) a combination of epoxide and hydroxyl molecules. Carbon
molecules are shown to be green, oxygen is blue, and hydrogen is pink.

Figure 5.3 Single and 4 layer graphene oxide structures containing 12% bonded oxygen species.
Carbon atoms are gray, oxygen is red, and hydrogen is white. Oxygen atoms protrude from the top
and bottom of the single layer, and the outermost layers of the 4 layer structures.

51
 After creating the models of graphene oxide and adding the polymer atoms as
described above. The models were densified in the same manner as in Chapter 4. All
models were condensed to very near the target density, and were then crosslinked.
The crosslinking procedure was the same as in previous works and is described in
detail in section 2.3 (page 12). All models were crosslinked to 75% crosslink density,
as results from the previous study for pristine graphite revealed little effect for
varying the amount of crosslinks. Fully equilibrated crosslinked structures for 1 and
4 layer graphene with the 3 different oxygen species are shown in Figure 5.4 for 12%
bonded oxygen. Before deforming the equilibrated models, the pressure was relived
using an NPT simulation causing an increase in the simulation box size for all 3
dimensions. In order to keep the density near the target, a final MD densification was
performed along the z-axis. Density values for each model can be seen in Table 5.1

Figure 5.4 Fully equilibrated Single and 4 layer graphene oxide / epoxy structures containing 12%
bonded oxygen species.

52
 Deformation simulations were then performed in both Tension and Shear, in
the exact manner as the previous chapter. Models were deformed to 5% strain,
providing tensile and shear modulus data, as well as Poisson’s ratio. Exact details are
provided in section 4.2 (page 34). All mechanical properties derived from the MD
models can be seen in Table 5.1 below. After deforming, results from the epoxide
only models were implemented into MAC/GMC software and compared to the
results for 75% crosslinked pristine GNP structures of 1 and 4 atomic thicknesses.

5.3 Results and Discussion

The interfacial polymer mass density profiles for the oxidized graphene show
less fluctuation near the graphene surface, as shown in Figure 5.5. This is likely due
to the fact that atoms are able to protrude between the strips of attached oxygen
species, making them appear more even. When examining the oxidized graphene
structures in Figure 5.2, each attached oxygen species causes a different amount of
distortion to the graphene surface. The amount of distortion (waviness) seen in each
of the models is related to the energy of the oxidized structure. As would be
expected, the epoxide only models show the least amount of distortion and the peaks
in mass density profile are easily visible for each layer. The hydroxyl models are
very wavy, causing the layers to be less distinguishable in the figure. This wavy
behavior is synonymous with buckling and proved to be a problem during
deformation, highlighted later.

53
Figure 5.5 Polymer mass density profile viewed along the Z-Axis for single layer a) 12% bonded
oxygen b) 24% bonded oxygen and 4 layer c) 12% bonded oxygen d) 24% bonded oxygen

Upon noticing the wavy behavior is the hydroxyl and combination models,
the per atom potential energy values were investigated. Figure 5.6 shows the energy
values for each atom within the graphene oxide structures for single and 4 layer
graphene oxide structures with varying oxygen species at 12% oxygen bonding. The
epoxide atoms seem to produce the highest energy values. This seems to go against
the referenced values from DFT, though those were calculated for the graphene
oxide as a whole rather than just the epoxide groups alone. In terms of bond energy,
the epoxide bond is a high energy bond, which is why it is a common reaction site.
Hydroxyl atoms cause the graphene to undulate, but are at a relatively low energy
state. The highest energy contributions from the combination graphene oxide models
are due to the epoxide molecules. The high energy bonding of the epoxide oxygen
molecules appears to contribute to the more planar structure. The undulated behavior
seen from the hydroxyl molecules results from the changing of the Carbon to Carbon
bond distances in the aromatic rings to which these molecules are bonded.

54
Figure 5.6 Visual representation of per-atom potential energy values (Kcal/mol-Angstrom) for single
and 4 layer graphite oxide structures for varying oxygen species with 12% oxygen bonding.

Mechanical properties are derived from stress vs. strain curves that look very
similar to those shown in Figure 4.5 from Chapter 4. All mechanical properties for
oxidized graphene models are shown in Table 5.1. Data from the 75% crosslinked
pristine graphene structures from Chapter 4 are also shown for comparison. When
allowing for Poisson’s contractions during tensile deformation, the wavy hydroxyl
group containing structures extended outward rather than contracting. This occurs
even after relieving the pressure along the x and y axes with a 2 nanosecond NPT
simulation. The Poisson’s values for which this phenomenon occurred simply have a
(+) value in the table, as the deformation of the oxidized graphene sheets in the axial
direction was coupled by extension transverse direction. This may be attributed to
the pattern in which the oxygen species were attached to the graphene surface. The
addition of oxygen species reduces the x and y axis (axial) tensile modulus values
when compared to that for pristine graphene structures. This effect is amplified by
55
adding more oxygen. The z axis (transverse) modulus values show an increase for
the oxygenated graphene structures when compared to pristine graphene, likely due
to stronger atomic interactions between the attached oxygen species and the
surrounding polymer atoms. Results for shear modulus along the xy plane are
decreased for the hydroxyl and combination models, though the epoxide models
produce slightly better results than that of pristine graphite. Unlike the results for
pristine graphite, shear in the xz direction provided positive values for all of the
models. The attached oxygen groups prevented the polymer atoms from easily
sliding past the graphite surface. The models utilizing only epoxide oxygen species
produced the best results for mechanical properties, including good results for
Poisson’s contractions and transverse isotropy, which would enable them to be
implemented into the MAC/GMC software.

Table 5.1 Mechanical properties derived from MD for oxidized graphene models along with those for
75% crosslinked pristine graphene MD models from chapter 4 for comparison

MAC/GMC results for oxidized graphene composed of only epoxide species

were very close to that for pristine graphite with 75% crosslinked epoxy, indicating
little influence from increases in mechanical properties shown in Table 5.1.
MAC/GMC micromechanics results are shown for up to 5% GNP concentration in
Table 5.2. As seen from the table, pristine graphite structures produce best results.

56
Table 5.2 MAC/GMC results for 75% crosslinked pristine and epoxide oxidized graphite structures of
1 and 4 atomic layers.

5.4 Conclusions

This chapter investigates the addition of oxygen species to graphene to

simulate graphene oxide in the presence of interfacial EPON 862-DETDA epoxy.
Even after an extensive relaxation process, some models did not accurately reflect
Poisson’s contractions during tensile deformation. Perhaps attaching the oxygen in a
different pattern would eliminate this problem. Interfacial characteristics agree with
those for pristine graphene, though the oscillations in polymer mass density profiles
are lessened due to the protruding hydroxyl molecules from the graphene surface.
Per-atom potential energy values are highest in the epoxide molecules, due to the
high energy of the epoxide bond. The addition of hydroxyl groups causes the
graphene to become wavy, though this may be due to the pattern of which they were
applied. Both elastic and shear modulus values are decreased with the addition of
hydroxyl oxygen to the graphene structures, and more oxygen appears to amplify this
effect. Epoxide only models are in best agreement with pristine graphene, followed
by models with a combination of epoxide and hydroxyl molecules. Attaching only

57
hydroxyl groups to the graphene surface produces the lowest values of tensile and
shear modulus. The addition of oxygen species provides increases in values for z-
axis tensile modulus and shear modulus in the x-z plane, though these increases have
little effect on the micromechanics results for a GNP doped composite.

CHAPTER 6: SCOPE FOR FUTURE WORKS

This chapter highlights the possibilities for future research in the realm of
epoxy / carbon composites. Suggestions for improvements to the research in this
report are also highlighted.

6.1 Improving Modeling Procedures

Throughout the research performed in this report, many possible changes to

the modeling procedure were brought to attention. The area for the most
improvement would be the densification equilibration and crosslinking procedures.
Crosslinking in this work was performed after the material was brought very near to
the target density of 1.17 g/cc. During the last research performed on the oxidized
graphene, all the models were brought near to the target density before crosslinking.
It was noticed that some models were able to achieve the target crosslink density of
75% much quicker than others. After investigating this, it was found that the models
that were further away from the target density were able to crosslink much quicker.
The reason for discrepancies in density from one model to another was that the
oxygen groups attached to the graphene surfaces would produce different amounts of
space between the graphene layers, along with the surrounding epoxy. More
specifically, models were densified based strictly on the number of layers of
graphene and the type of oxygen species attached to the surface was disregarded.
58
This produced a different density in each model prior to crosslinking. When these
materials crosslink, there is a small amount of volume shrinkage. Based on this, it
could be proposed that models be brought very near the target and then crosslinked.
The concept of volume shrinkage could be considered and employed during the
densification/crosslinking procedure. With the LAMMPS software always
improving, it is now possible to crosslink and densify at the same time all within the
MD simulation. Performing these procedures statically outside of the actual MD
simulations is very time consuming and if eliminated could vastly speed up the
process. Not to mention, it is more representative of a realistic cross-linking and
densification method. It is also worth mentioning that these materials are often cured
at very high temperatures, which could easily be accounted for with MD. In
summary, the crosslinking and densification process could be employed at the same
time with temperature considerations specific to the material being studied. The
overall effect on material characteristics would likely be minimal, but the suggested
method would be less subject to question.

Another improvement to the models discussed in this report would be the

consideration of atomic charges within the model. Atomic Charges would influence
the crosslinking mechanics during the simulations, along with the behavior of the
oxygen species attached to the oxidized graphene surfaces. Again, the overall
influence on the end product would likely be minimal.

6.2 Modeling Functionalized Graphene and Surface Treatments

The existence of oxygen molecules on the surface of graphite structures

enables for the bonding of larger molecules, which could then bond to the
surrounding epoxy/hardener molecules. This process is known as functionalization,
and has been seen to produce significant improvements in composite properties.65-67
Molecular Dynamics models have been developed for a variety of functionalization

59
methods.63,64,68-70 To fully understand the interactions between polymer matrix and
carbon reinforcement, models should be developed for every possible scenario that
can occur in nature. As the interfacial mechanisms become more understood,
models can be developed that incorporate all of the interactions that may occur (non-
bonded Van der Waals forces, charges from oxygen and nitrogen species, and
bonding between carbon surface and surrounding polymer matrix). As the models
become more accurate (realistic), they will become a more reliable tool for
developing materials for very specific applications.

The use of MD could also be used to model various surface treatments that
are performed during fabrication of polymer/carbon composites. Experimental
research has revealed certain treatments to be beneficial to overall composite
mechanical properties. These include: ultrasonic treatment71, plasma oxidation72,
electron beam irradiation73, DQGȖ-ray irradiation74. The reason that these treatments
work is because they alter the molecular behavior of the materials, something that
can only be captured by MD.

6.3 Modeling Aging in Composites

Epoxy composites can often be subject to environments that promote both

physical and chemical aging. Physical aging occurs when materials are exposed to
temperatures below the glass transition for extended periods of time. This can cause
a reduction in creep for composite materials.75 Chemical aging involves
photochemical, oxidative, or thermal degradation which occurs due to formation or
breaking of covalent bonds within a material. The most common form of chemical
aging for composites is hygrothermal, which occurs in moist environments at high
temperatures.76 Hygrothermal aging severs chemical bonds through hydrolysis and
can result in loss of adhesion between epoxy and fiber.

60
 Molecular Dynamics has been used to model physical aging of the EPON
862/DETDA polymer system20, though it has never been investigated for a
composite material. MD could be particularly useful in investigating the
hygrothermal degradation process for an epoxy composite. It would be relatively
easy to develop a composite model and then insert water molecules and increase the
simulation temperature. This could be done for a variety of fiber/matrix systems.

The possibilities for research in the realm of MD are endless. Models will
continue to become larger and more accurate as computational power increases and
experimental techniques continue to reveal new information.

61
 REFERENCES
1. Drzal, L.T. and M. Madhukar, Fiber Matrix Adhesion and Its Relationship to
Composite Mechanical-Properties. Journal of Materials Science, 1993. 28(3): p.
569-610.
2. Wagner, D.H. and R.A. Vaia, Nanocomposites: Issues at the Interface. Materials
Today, 2004. 7(11): p. 38-42.
3. Potter, W. G., Epoxide Resins. In Springer-Verlag: New York, 1970.
4. Bascom, W.D. and Drzal, L.T., The Surface Properties of Carbon Fibers and Their
Adhesion to Organic Polyers. NASA Contractor Report 4084, July 1987.
5. Drzal, L., The interphase in epoxy composites, in Epoxy Resins and Composites II,
K. Dušek, Editor. 1986, Springer Berlin Heidelberg. p. 1-32.
6. Gojny, F.H., et al., Carbon nanotube-reinforced epoxy-composites: enhanced
stiffness and fracture toughness at low nanotube content. Composites Science and
Technology, 2004. 64(15): p. 2363-2371.
7. Gojny, F.H., et al., Influence of different carbon nanotubes on the mechanical
properties of epoxy matrix composites – A comparative study. Composites Science
and Technology, 2005. 65(15–16): p. 2300-2313.
8. RamanathanT, et al., Functionalized graphene sheets for polymer nanocomposites.
Nat Nano, 2008. 3(6): p. 327-331.
9. Rafiee, M.A., et al., Enhanced Mechanical Properties of Nanocomposites at Low
Graphene Content. ACS Nano, 2009. 3(12): p. 3884-3890.
10. King, J.A., et al., Mechanical properties of graphene nanoplatelet/epoxy composites.
Journal of Applied Polymer Science, 2013. 128(6): p. 4217-4223.
11. Mäder, E. and E. Pisanova, Interfacial design in fiber reinforced polymers.
Macromolecular Symposia, 2001. 163(1): p. 189-212.
12. LAMMPS website: http://lammps.sandia.gov
13. Stevens, M. J. Macromolecules 2001, 34, (8), 2710-2718.
14. Tsige, M.; Stevens, M. J. Macromolecules 2004, 37, (2), 630-637.
15. Carmesin, I.; Kremer, K. Macromolecules 1988, 21, (9), 2819-2823.
16. Jo, W. H.; Ko, M. B. Macromolecules 1994, 27, (26), 7815-7824.
17. Jo, W. H.; Ko, M. B. Macromolecules 1993, 26, (20), 5473-5478.

62
18. Fan, H.B. and M.M.F. Yuen, Material properties of the cross-linked epoxy resin
compound predicted by molecular dynamics simulation. Polymer, 2007. 48(7): p.
2174-2178.
19. Varshney, V., S. Patnaik, A. Roy, and B. Farmer, A Molecular Dynamics Study of
Epoxy Based Networks: Cross-linking Procedure and Prediction of Molecular and
Material Properties. Macromolecules, 2008. 41(18): p. 6837-6842
20. Bandyopadhyay, A., P.K. Valavala, T.C. Clancy, K.E. Wise, and G.M. Odegard,
Molecular modeling of crosslinked epoxy polymers: The effect of crosslink density
on thermomechanical properties. Polymer, 2011. 52(11): p. 2445-2452
21. Li, C.Y. and A. Strachan, Molecular dynamics predictions of thermal and
mechanical properties of thermoset polymer EPON862/DETDA. Polymer, 2011.
52(13): p. 2920-2928
22. Frankland, S.J.V., A. Caglar, D.W. Brenner, and M. Griebel, Molecular Simulation
of the Influence of Chemical Cross-Links on the Shear Strength of Carbon
Nanotube-Polymer Interfaces. Journal of Physical Chemistry B, 2002. 106(12): p.
3046-3048
23. Odegard, G.M., T.S. Gates, K.E. Wise, C. Park, and E.J. Siochi, Constitutive
modeling of nanotube-reinforced polymer composites. Composites Science and
Technology, 2003. 63(11): p. 1671-1687.
24. Frankland, S.J.V., V.M. Harik, G.M. Odegard, D.W. Brenner, and T.S. Gates, The
Stress-Strain Behavior of Polymer-Nanotube Composites from Molecular Dynamics
Simulation. Composites Science and Technology, 2003. 63(11): p. 1655-1661.
25. Clancy, T.C. and T.S. Gates, Modeling of interfacial modification effects on thermal
conductivity of carbon nanotube composites. Polymer, 2006. 47(16): p. 5990-5996.
26. Zhu, R., E. Pan, reinforced Epon 862 composites. Materials Science and
Engineering a-Structural Materials Properties and A.K. Roy, Molecular dynamics
study of the stress-strain behavior of carbon-nanotube Microstructure and
Processing, 2007. 447(1-2): p. 51-57.
27. Nouranian, S., C. Jang, T.E. Lacy, S.R. Gwaltney, H. Toghiani, and C.U. Pittman,
Molecular dynamics simulations of vinyl ester resin monomer interactions with a
pristine vapor-grown carbon nanofiber and their implications for composite
interphase formation. Carbon, 2011. 49(10): p. 3219-3232.

63
28. Gou, J., et al., Computational and experimental study of interfacial bonding of
single-walled nanotube reinforced composites. Computational Materials Science,
2004. 31(3–4): p. 225-236.
29. Odegard, G.M., T.C. Clancy, and T.S. Gates, Modeling of the mechanical properties
of nanoparticle/polymer composites. Polymer, 2005. 46(2): p. 553-562.
30. Yu, S., S. Yang, and M. Cho, Multi-scale modeling of cross-linked epoxy
nanocomposites. Polymer, 2009. 50(3): p. 945-952.
31. Jang, C., et al., Interfacial shear strength of cured vinyl ester resin-graphite
nanoplatelet from molecular dynamics simulations. Polymer, 2013. 54(13): p. 3282-
3289
32. Mansfield, K.F. and D.N. Theodorou, Atomistic Simulation of a Glassy Polymer
Graphite Interface. Macromolecules, 1991. 24(15): p. 4295-4309
33. Stevens, M.J., Interfacial Fracture between Highly Cross-Linked Polymer Networks
and a Solid Surface: Effect of Interfacial Bond Density. Macromolecules, 2001. 34:
p. 2710-2718.
34. Alkhateb, H., A. Al-Ostaz, and A. H. D. Cheng. "Molecular dynamics simulations of
graphite-vinyl ester nanocomposites and their constituents." Carbon Lett 11.4
(2010): 316-324.
35. Hadden, C.M., et al., Molecular modeling of EPON-862/graphite composites:
Interfacial characteristics for multiple crosslink densities. Composites Science and
Technology, 2013. 76(0): p. 92-99.
36. Gao, J.-S., S.-C. Shiu, and J.-L. Tsai, Mechanical properties of polymer near
graphite sheet. Journal of Composite Materials, 2013. 47(4): p. 449-458.
37. Li, C., et al., Atomistic simulations on multilayer graphene reinforced epoxy
composites. Composites Part A: Applied Science and Manufacturing, 2012. 43(8): p.
1293-1300.
38. Mittal, R., et al. Molecular modeling for calculation of mechanical properties of
EPON862 / SWCNTs composites. in Nanoscience, Engineering and Technology
(ICONSET), 2011 International Conference on. 2011.
39. Han, Y. and J. Elliott, Molecular dynamics simulations of the elastic properties of
polymer/carbon nanotube composites. Computational Materials Science, 2007.
39(2): p. 315-323.

64
40. Ionita, M., Multiscale molecular modeling of SWCNTs/epoxy resin composites
mechanical behaviour. Composites Part B: Engineering, 2012. 43(8): p. 3491-3496.

41. Duffy, E. M.; Kowalczyk, P. J.; Jorgensen, W. L. Journal of the American Chemical
Society 1993, 115, (20), 9271-9275.
42. Jorgensen, W. L.; Maxwell, D. S.; TiradoRives, J. Journal of the American
Chemical Society 1996, 118, (45), 11225-11236.
43. Weiner, S. J.; Kollman, P. A.; Case, D. A.; Singh, U. C.; Ghio, C.; Alagona, G.;
Profeta, S.; Weiner, P. Journal of the American Chemical Society 1984, 106, (3),
765-784.
44. Watkins, E. K.; Jorgensen, W. L. Journal of Physical Chemistry A 2001, 105, (16),
4118-4125.
45. Plimpton, S., Fast Parallel Algorithms for Short-Range Molecular-Dynamics
Journal of Computational Physics, 1995. 117(1): p. 1-19.
46. Hoover, W.G., Canonical Dynamics - Equilibrium Phase-Space Distributions.
Physical Review A, 1985. 31(3): 1695-1697.
47. Yarovsky, I. and E. Evans, Computer simulation of structure and properties of
crosslinked polymers: application to epoxy resins. Polymer, 2002. 43(3): 963-969.
48. Littell, J.D., C.R. Ruggeri, R.K. Goldberg, G.D. Roberts, W.A. Arnold, and W.K.
Binienda, Measurement of epoxy resin tension, compression, and shear stress-strain
curves over a wide range of strain rates using small test specimens. Journal of
Aerospace Engineering, 2008. 21(3): p. 162-173.
49. Kuilla, T., et al., Recent advances in graphene based polymer composites. Progress
in Polymer Science, 2010. 35(11): p. 1350-1375.
50. Tang, L.-C., et al., The effect of graphene dispersion on the mechanical properties of
graphene/epoxy composites. Carbon, 2013. 60(0): p. 16-27
51. Aboudi, J.; Arnold, S. M.; Bednarcyk, B. A. (2013): Micromechanics of Composite
Materials: A Generalized Multiscale Analysis Approach. Elsevier, Inc.
52. Bednarcyk, B.A., J. Aboudi, and S.M. Arnold, The effect of general statistical fiber
misalignment on predicted damage initiation in composites. Composites Part B:
Engineering, 2014. 66(0): p. 97-108.
53. Mechanics of Materials, 1992. 14(2): p. 127-139.Paley, M. and J. Aboudi,
Micromechanical analysis of composites by the generalized cells model.

65
54. HexTow ® Continuous Carbon Fiber Product Literature; Hexcel: Stamford, CT
(2011)
55. Eda, G., et al., Partially oxidized graphene as a precursor to graphene. Journal of
Materials Chemistry, 2011. 21(30): p. 11217-11223.
56. Szabó, T., et al., Evolution of Surface Functional Groups in a Series of
Progressively Oxidized Graphite Oxides. Chemistry of Materials, 2006. 18(11): p.
2740-2749
57. Dreyer, D.R., et al., The chemistry of graphene oxide. Chemical Society Reviews,
2010. 39(1): p. 228-240.

58. Allington, R.D., et al., A model of the surface of oxidatively treated carbon fibre
based on calculations of adsorption interactions with small molecules. Composites
Part A: Applied Science and Manufacturing, 1998. 29(9–10): p. 1283-1290.
59. Allington, R.D., et al., Developing improved models of oxidatively treated carbon
fibre surfaces, using molecular simulation. Composites Part A: Applied Science and
Manufacturing, 2004. 35(10): p. 1161-1173.
60. Kudin, K.N., et al., Raman Spectra of Graphite Oxide and Functionalized Graphene
Sheets. Nano Letters, 2007. 8(1): p. 36-41.
61. Boukhvalov, D.W. and M.I. Katsnelson, Modeling of Graphite Oxide. Journal of the
American Chemical Society, 2008. 130(32): p. 10697-10701.
62. Sadiq Awan, I., et al., Developing an approach to calculate carbon fiber surface
energy using molecular simulation and its application to real carbon fibers. Journal
of Composite Materials, 2011.
63. Xiaoqun, W., et al., Molecular dynamics simulations of interfacial adhesion between
carbon fibers and various epoxies/hardeners and its calorimetric validation. Journal
of Composite Materials, 2013. 47(8): p. 1011-1017.
64. Rahman, R. and A. Haque, Molecular modeling of crosslinked graphene–epoxy
nanocomposites for characterization of elastic constants and interfacial properties.
Composites Part B: Engineering, 2013. 54(0): p. 353-364.

65. Ganguli, S., A.K. Roy, and D.P. Anderson, Improved thermal conductivity for
chemically functionalized exfoliated graphite/epoxy composites. Carbon, 2008.
46(5): p. 806-817.

66
66. Miller, S.G., et al., Characterization of epoxy functionalized graphite nanoparticles
and the physical properties of epoxy matrix nanocomposites. Composites Science
and Technology, 2010. 70(7): p. 1120-1125.

67. Mäder, E. and E. Pisanova, Interfacial design in fiber reinforced polymers.

Macromolecular Symposia, 2001. 163(1): p. 189-212

68. Attwood, D. and P.I. Marshall, Atomistic modelling of the adsorption of epoxy and
amine molecules on the surface of carbon fibres. Composites Part A: Applied
Science and Manufacturing, 1996. 27(9): p. 775-779.
69. Sanchez, F. and L. Zhang, Molecular dynamics modeling of the interface between
surface functionalized graphitic structures and calcium–silicate–hydrate:
Interaction energies, structure, and dynamics. Journal of Colloid and Interface
Science, 2008. 323(2): p. 349-358.
70. Aparicio, S. and M. Atilhan, Choline-Based Ionic Liquids on Graphite Surfaces and
Carbon Nanotubes Solvation: A Molecular Dynamics Study. The Journal of Physical
Chemistry C, 2012. 116(22): p. 12055-12065.

71. Huang, Y.D., et al., Influence of ultrasonic treatment on the characteristics of epoxy
resin and the interfacial property of its carbon fiber composites. Composites
Science and Technology, 2002. 62(16): p. 2153-2159.
72. Montes-Morán, M.A., et al., Effects of plasma oxidation on the surface and
interfacial properties of carbon fibres/polycarbonate composites. Carbon, 2001.
39(7): p. 1057-1068.

73. Zhang, Z., et al., The effect of carbon-fiber surface properties on the electron-beam
curing of epoxy-resin composites. Composites Science and Technology, 2002. 62(3):
p. 331-337.
74. Xu, Z., et al., (IIHFWRIȖ-ray irradiation grafting on the carbon fibers and interfacial
adhesion of epoxy composites. Composites Science and Technology, 2007. 67(15–
16): p. 3261-3270.
75. Sullivan, J.L., Creep and physical aging of composites. Composites Science and
Technology, 1990. 39(3): p. 207-232
76. Ray, B.C., Temperature effect during humid ageing on interfaces of glass and
carbon fibers reinforced epoxy composites. Journal of Colloid and Interface Science,
2006. 298(1): p. 111-117.
67
APPENDIX A
LAMMPS scripts for MD simulations, MM simulations, and deformation
simulations

A.1. Molecular Energy Minimization Script

#Energy minimisation input script

units real
dimension 3
boundary p p p
atom_style molecular
neighbor 2.0 bin
neigh_modify one 25000 page 250000
echo screen

#OPLS potentials
bond_style harmonic
angle_style harmonic
dihedral_style opls
pair_style lj/cut/opt 10.0
pair_modify mix arithmetic
read_data filename.xyz

#GROUPS AND FIXES (FROM MODELS IN CHAPTER 3)

group plug id > 21120
group graphite type 14
fix 1 plug setforce 0.0 0.0 0.0
fix 2 graphite setforce 0.0 0.0 0.0
neigh_modify exclude group plug plug
neigh_modify exclude group graphite graphite
dump 1 all xyz 500 dump.filename.xyz
dump 2 all atom 500 dump.filename.lammpstrj
dump_modify 2 scale no

#minimisation
min_style cg #conjugate gradient
minimize 0.0 0.0 10000 100000
write_restart min.filename.restart

68
A.2. Molecular Dynamics NVT simulation script for equilibration

# NVT simulation for equilibration

units real
dimension 3
boundary p p p
atom_style molecular
neighbor 5.0 bin
neigh_modify every 1 delay 10 one 4000 page 250000
echo screen

#----------------OPLS potentials--------------------
bond_style harmonic
angle_style harmonic
dihedral_style opls
pair_style lj/cut/opt 10.0
pair_modify mix arithmetic
read_restart min.filename.restart

#---------Define Groups for specific output---------

group graphite type 14 16 17 18
group benzene type 1 8
#group xlink type 16 17 18 19 20

#------Commands for freezing atoms---------------

#neigh_modify exclude group plug plug
#neigh_modify exclude group graphite graphite

#--------------initial velocities------------------------------
velocity all create 300.0 3094152
compute 1 all temp
compute 2 all pe

#-----per-atom potential energy and virial contributions-------

compute spa all stress/atom
compute sta all reduce sum c_spa[1] c_spa[2] c_spa[3] c_spa[4] c_spa[5]
c_spa[6]
variable spaxx atom c_spa[1]/vol*(1/9869.23267) #GPa
variable spayy atom c_spa[2]/vol*(1/9869.23267)
variable spazz atom c_spa[3]/vol*(1/9869.23267)
variable spaxy atom c_spa[4]/vol*(1/9869.23267)
variable spaxz atom c_spa[5]/vol*(1/9869.23267)
variable spayz atom c_spa[6]/vol*(1/9869.23267)
variable staxx equal c_sta[1]/vol*(1/9869.23267)

69
variable stayy equal c_sta[2]/vol*(1/9869.23267)
variable stazz equal c_sta[3]/vol*(1/9869.23267)
variable staxy equal c_sta[4]/vol*(1/9869.23267)
variable staxz equal c_sta[5]/vol*(1/9869.23267)
variable stayz equal c_sta[6]/vol*(1/9869.23267)
variable spaVM atom sqrt((v_spaxx-v_spayy)^2+(v_spayy-
v_spazz)^2+(v_spaxx-v_spazz)^2+6*(v_spaxy^2+v_spaxz^2+v_spayz^2)) # von
Mises stress
compute ppa all pe/atom
compute pta all reduce ave c_ppa
variable ppa atom c_ppa
#compute ppavdw all pe/atom pair
variable ppavdw atom c_ppavdw
#compute ppabnd all pe/atom bond
#compute ppaang all pe/atom angle
#compute ppadih all pe/atom dihedral

#-------fixes for calculating various denity, stress, and energy contributions------

fix 1 all nvt temp 300.0 300.0 10.0
#--------density calculations along z-axis----------------
fix 2 all ave/spatial 500 50 500000 z lower 0.1 density/mass units box ave one file
density1.504.c75.eq.txt
#--------densification command---------------------------
fix 3 all deform 1 z final -16.15 16.15 units box
#--------specific outputs - alter as desired---------------
#fix 5 benzene ave/spatial 500 50 25000 z lower 0.2 density/mass units box ave one
file density.filename.benzene.txt
#fix 6 all ave/spatial 1 1 1 z lower 0.1 c_ppa c_ppavdw c_ppabnd c_ppaang
c_ppadih units box ave one file energy.filename.txt
#fix 6 all ave/spatial 500 50 25000 z lower 0.1 c_ppa c_ppavdw c_ppabnd c_ppaang
c_ppadih units box ave one file energy.filename.txt
#fix 7 all ave/spatial 1 1 1 z lower 0.1 v_spaxx v_spayy v_spazz v_spaxz v_spaxy
v_spayz v_spalen units box ave one file stress.filename.txt
fix 4 all ave/spatial 500 50 50000 z lower 0.1 v_spaxx v_spayy v_spazz v_spaxz
v_spaxy v_spayz v_spalen units box ave one file stress.filename.txt
#fix 8 benzene ave/spatial 500 50 25000 z lower 0.1 v_spaxx v_spayy v_spazz
v_spaxz v_spaxy v_spayz v_spalen units box ave one file stress.filename.benzene.txt

#----OVITO stress and energy dump files for visualization---

#dump normal all custom 500000 stress.filename.lammpstrj id type x y z v_spaxx
v_spayy v_spazz v_spaVM
#dump_modify normal sort id flush yes
#dump shear all custom 5000 shear.filename.lammpstrj id type x y z v_spaxy
v_spaxz v_spayz
#dump_modify shear sort id flush yes

70
dump pe_st graphite custom 50000 Graphite.energy.lammpstrj id type x y z v_ppa
dump_modify pe_st sort id flush yes

#-------------------run settings-----------------------------
timestep 0.2 #femtosecons
thermo 500 # 0.1 picoseconds
thermo_style custom step temp press pxx pyy pzz pe ke etotal evdwl ebond eangle
edihed emol vol enthalpy
#Old Methods
dump 1 graphite atom 50000 dump.Graphite.lammpstrj
dump_modify 1 scale no
dump 2 all xyz 500000 dump.filename.xyz #xyz files are necessary for use of
crosslinking program
run 2000000 #0.4 nanosecond
write_restart MD.nvt1.504.c75.eq.restart

71
A.3. Molecular Dynamics NPT simulation script for Tensile Deformation
# NPT Simulation for Tensile Deformation

#-----------------Initialization----------------------
units real
dimension 3
boundary p p p
atom_style molecular
neighbor 5.0 bin
neigh_modify every 1 delay 10 one 4000 page 250000
echo screen

#----------------OPLS potentials--------------------
bond_style harmonic
angle_style harmonic
dihedral_style opls
pair_style lj/cut 10.0
pair_modify mix arithmetic
read_restart min.filename.Exx.restart
group graphite type 14 16 17 18
group polymer type 1 2 3 4 5 6 7 8 9 10 11 12 13 15 19 20 21 22
group xlink type 16 17 18 19 20

#--------------initial velocities------------------------------
compute 1 all temp
compute 2 all pe
velocity all create 300.0 300430

#-----per-atom potential energy and virial contributions------

variable Time equal step*dt/1000 #time in picoseconds
variable mass_g equal mass(all)/6.02e23
variable vol_cm3 equal vol*1.0e-24
variable massdensity equal v_mass_g/v_vol_cm3
variable lx equal lx
variable ly equal ly
variable lz equal lz
variable lx0 equal ${lx}
variable ly0 equal ${ly}
variable lz0 equal ${lz}
variable esxx equal (v_lx-v_lx0)/v_lx0
variable esyy equal (v_ly-v_ly0)/v_ly0
variable eszz equal (v_lz-v_lz0)/v_lz0
variable tsxx equal ln(v_lx/v_lx0)
variable tsyy equal ln(v_ly/v_ly0)
variable tszz equal ln(v_lz/v_lz0)

72
variable pxx equal pxx
variable pyy equal pyy
variable pzz equal pzz
variable tsigxx equal -(v_pxx)*0.101325 #true stress xx
variable tsigyy equal -(v_pyy)*0.101325 #true stress yy
variable tsigzz equal -(v_pzz)*0.101325 #true stress zz
compute spa all stress/atom
compute sta all reduce sum c_spa[1] c_spa[2] c_spa[3] c_spa[4] c_spa[5]
c_spa[6]
variable spaxx atom c_spa[1]/vol/9.86923
variable spayy atom c_spa[2]/vol/9.86923
variable spazz atom c_spa[3]/vol/9.86923
variable spaxy atom c_spa[4]/vol/9.86923
variable spaxz atom c_spa[5]/vol/9.86923
variable spayz atom c_spa[6]/vol/9.86923
variable staxx equal c_sta[1]/vol/9.86923
variable stayy equal c_sta[2]/vol/9.86923
variable stazz equal c_sta[3]/vol/9.86923
variable staxy equal c_sta[4]/vol/9.86923
variable staxz equal c_sta[5]/vol/9.86923
variable stayz equal c_sta[6]/vol/9.86923
variable spalen atom
sqrt((v_staxx*v_staxx)+(v_stayy*v_stayy)+(v_stazz*v_stazz))

#----------------Graphite Computes----------------------
compute Gspa graphite stress/atom
compute Gsta graphite reduce sum c_Gspa[1] c_Gspa[2] c_Gspa[3] c_Gspa[4]
c_Gspa[5] c_Gspa[6]
variable Gspaxx atom c_Gspa[1]/vol/9.86923
variable Gspayy atom c_Gspa[2]/vol/9.86923
variable Gspazz atom c_Gspa[3]/vol/9.86923
variable Gstaxx equal c_Gsta[1]/vol/9.86923
variable Gstayy equal c_Gsta[2]/vol/9.86923
variable Gstazz equal c_Gsta[3]/vol/9.86923

#-------------------polymer computes--------------------------
compute Pspa polymer stress/atom
compute Psta polymer reduce sum c_Pspa[1] c_Pspa[2] c_Pspa[3] c_Pspa[4]
c_Pspa[5] c_Pspa[6]
variable Pspaxx atom c_Pspa[1]/vol/9.86923
variable Pspayy atom c_Pspa[2]/vol/9.86923
variable Pspazz atom c_Pspa[3]/vol/9.86923
variable Pstaxx equal c_Psta[1]/vol/9.86923
variable Pstayy equal c_Psta[2]/vol/9.86923

73
variable Pstazz equal c_Psta[3]/vol/9.86923

#-------------------Energy Contributions-----------------------
compute kpa all ke/atom #per-atom kinetic energy
compute ppa all pe/atom #per-atom potential energy
compute pta all reduce ave c_ppa #total potential energy
variable kpa equal c_kpa
variable ppa atom c_ppa
variable pta equal c_pta

compute Gppa graphite pe/atom #per-atom potential energy

compute Gpta graphite reduce ave c_Gppa#total PE
variable Gppa equal c_Gppa
variable Gpta equal c_Gpta

compute Pppa polymer pe/atom #per-atom potential energy

compute Ppta polymer reduce ave c_Pppa #total pot energy
variable Pppa equal c_Pppa
variable Ppta equal c_Ppta
#thermo_style custom v_Time temp v_massdensity press pxx pyy pzz pe v_tsxx
v_tsyy v_tszz v_tsigxx v_tsigyy v_tsigzz

#compute ppavdw all pe/atom pair

variable ppavdw atom c_ppavdw
#compute ppabnd all pe/atom bond
#compute ppaang all pe/atom angle
#compute ppadih all pe/atom dihedral

#------------MAIN DEFORMATION COMMAND USE NPT FOR TENSILE-----

fix 1 all npt temp 300.0 300.0 100.0 y 1.0 1.0 1000.0 z 1.0 1.0 1000.0 couple none
fix 2 all deform 1 x scale 1.05 remap none units box #DEFORM TO 5% STRAIN
#-TENSILE DATA OUTPUT COMMAND FOR MAKING STRESS STRAIN
CURVES--
fix 7 all ave/time 1 100 1000 v_esxx v_esyy v_eszz v_staxx v_stayy v_stazz file
tensile.filename.txt

#-----OVITO stress and energy dump files for visualization------

dump normal all custom 500000 Exx.stress.output.lammpstrj id type x y z v_spaxx
v_spayy v_spazz
dump_modify normal sort id flush yes
#dump shear all custom 5000 shear.filename.lammpstrj id type x y z v_spaxy
v_spaxz v_spayz
#dump_modify shear sort id flush yes
dump pe_st all custom 500000 Exx.energy.lammpstrj id type x y z v_ppa

74
dump_modify pe_st sort id flush yes
#-------------------run settings-----------------------------
thermo_style custom step v_Time temp vol pzz pxx pyy ke pe etotal v_pta v_Gpta
v_Ppta
timestep 0.2 #femtosecons
thermo 1000 # 0.2 picoseconds
#thermo_style custom step temp press pxx pyy pzz pe ke etotal evdwl ebond eangle
edihed emol vol enthalpy

#------------dump files for visualization---------

dump 1 all atom 500000 dump.nvt1.504.c75.Exx.lammpstrj
dump_modify 1 scale no
dump 2 all xyz 1000000 dump.nvt1.504.c75.Exx.xyz
run 5000000 # 1 nanosecond
write_restart deformed1.504.c75.Exx.restart

75
A.4. Molecular Dynamics NVT simulation script for Shear Deformation

# NVT simulation for Shear Deformation

#-----------------Initialization----------------------
Same as A.3
#----------------AMBER potentials--------------------
Same as A.3
#--------------initial velocities------------------------------
Same as A.3
#-----per-atom potential energy and virial contributions------
Same as A.3

#------MAIN DEFORMATION COMMAND USE NVT FOR SHEAR---------------

fix 1 all nvt temp 300.0 300.0 10.0
fix 5 all deform 1 xy final 5.0 remap none units box #DEFORM 5 ANGSTROMS
~5% OF X AND Y BOX SIZE ~100Angstroms
#---SHEAR DATA OUTPUT COMMAND FOR STRESS STRAIN CURVES---
fix 7 all ave/time 1 100 1000 v_xy v_yz v_xz v_staxy v_stayz v_staxz file
Gxy.shear.filename.txt

#------OVITO stress and energy dump files---------------

dump normal all custom 500000 Gxy.stress.lammpstrj id type x y z v_spaxx v_spayy
v_spazz v_spalen
dump_modify normal sort id flush yes
dump shear all custom 500000 Gxy.shear.lammpstrj id type x y z v_spaxy v_spaxz
v_spayz
dump_modify shear sort id flush yes
dump pe_st all custom 500000 Gxy.energy.lammpstrj id type x y z v_ppa
dump_modify pe_st sort id flush yes

#-------------------run settings-----------------------------
thermo_style custom step v_Time temp vol pzz pxx pyy ke pe etotal v_pta v_Gpta
v_Ppta
timestep 0.2 #femtosecons
thermo 1000 # 0.2 picoseconds
#thermo_style custom step temp press pxx pyy pzz pe ke etotal evdwl ebond eangle
edihed emol vol enthalpy

#-------------------Dump commands for visualization------------

dump 1 all atom 500000 dump.Gxy.filename.lammpstrj
dump_modify 1 scale no
dump 2 all xyz 1000000 dump.Gxy.filename.xyz
run 5000000 # 1 nanosecond
write_restart Gxy.filename.restart

76
A.5. MD OPLS UA Parameters from Coordinates File

#DATA FOR UNITED ATOM SYSTEM WITH 1 HYDROXYL AND EPOXIDE

FUNCTIONALIZED GRAPHITE SHEET

25857 atoms
23 atom types

29558 bonds
27 bond types

48419 angles
40 angle types

77086 dihedrals
56 dihedral types

-6.0599999999999998e-01 1.0053800000000000e+02 xlo xhi

-7.0999999999999996e-01 1.0356999999999999e+02 ylo yhi
-1.6199999999999999e+01 1.6199999999999999e+01 zlo zhi

Masses
atomid mass
1 12.011 #benzene carbon
2 1.008 #benzene hydrogen
3 14.027 #CH2 bonded to ring
4 14.027 #CH2 bonded to O
5 15.9994 #oxygen before x-link
6 13.017 #CH ether
7 15.9994 #oxygen bonded to ring
8 12.011 #carbon bonded to O
9 15.035 #CH3
10 14.007 #non-xlink nitrogen
11 1.008 #hydrogen x-link
12 12.0 #carbon from CH3
13 1.008 #hydrogen from CH3
14 12.0 #graphene carbon
15 14.027 #x-link CH2
16 15.9994 #functional hydroxyl O
17 1.008 #functional hydroxyl H
18 15.9994 #functional epoxide O
19 14.007 #nitrogen 1-xlink
20 15.9994 #crosslinked O
21 14.027 #crosslinked CH2
22 1.008 #crosslinked H

77
23 14.007 #nitrogen 2-xlinks

Bond Coeffs
bondtype K r0
16 320.0 1.425 #CH2 bonded to O - oxygen bonded to ring
10 317.0 1.510 #benzene carbon - CH2 bonded to ring
23 434.0 1.010 #non-xlink nitrogen - hydrogen x-link
21 320.0 1.425 #oxygen before x-link - CH ether
2 260.0 1.526 #CH ether - crosslinked CH2
18 434.0 1.010 #hydrogen x-link - nitrogen 1-xlink
8 367.0 1.080 #benzene carbon - benzene hydrogen
6 340.0 1.090 #carbon from CH3 - hydrogen from CH3
26 450.0 1.364 #graphene carbon - functional epoxide O
25 469.0 1.400 #graphene carbon - graphene carbon
22 481.0 1.340 #benzene carbon - non-xlink nitrogen
20 260.0 1.526 #CH ether - x-link CH2
7 469.0 1.400 #benzene carbon - benzene carbon
17 337.0 1.463 #nitrogen 1-xlink - crosslinked CH2
24 553.0 0.945 #functional hydroxyl O - functional hydroxyl H
5 553.0 0.945 #crosslinked O - crosslinked H
9 481.0 1.340 #benzene carbon - nitrogen 2-xlinks
4 337.0 1.463 #crosslinked CH2 - nitrogen 2-xlinks
11 481.0 1.340 #benzene carbon - nitrogen 1-xlink
27 450.0 1.364 #graphene carbon - functional hydroxyl O
15 260.0 1.526 #CH2 bonded to O - CH ether
12 469.0 1.400 #benzene carbon - carbon from CH3
3 386.0 1.425 #CH ether - crosslinked O
13 469.0 1.400 #benzene carbon - carbon bonded to O
1 260.0 1.526 #CH2 bonded to ring - CH3
19 320.0 1.425 #oxygen before x-link - x-link CH2
14 450.0 1.364 #oxygen bonded to ring - carbon bonded to O

Angle Coeffs
angletype K theta0
17 80.0 109.5 #CH2 bonded to O - CH ether - crosslinked O
37 35.0 118.4 #hydrogen x-link - nitrogen 1-xlink - crosslinked CH2
14 63.0 112.4 #benzene carbon - CH2 bonded to ring - benzene carbon
26 80.0 109.5 #crosslinked O - CH ether - crosslinked CH2
15 35.0 109.5 #benzene carbon - carbon from CH3 - hydrogen from CH3
38 43.6 106.4 #hydrogen x-link - non-xlink nitrogen - hydrogen x-link
20 80.0 109.5 #CH ether - CH2 bonded to O - oxygen bonded to ring
12 63.0 114.0 #benzene carbon - CH2 bonded to ring - CH3
1 85.0 120.0 #graphene carbon - graphene carbon - graphene carbon
32 35.0 111.0 #benzene carbon - nitrogen 1-xlink - hydrogen x-link
25 85.0 120.0 #benzene carbon - benzene carbon - carbon bonded to O

78
8 35.0 120.0 #benzene hydrogen - benzene carbon - carbon bonded to O
29 80.0 111.2 #CH ether - crosslinked CH2 - nitrogen 2-xlinks
23 63.0 111.1 #CH2 bonded to O - CH ether - x-link CH2
40 70.0 120.0 #benzene carbon - benzene carbon - non-xlink nitrogen
5 100.0 111.8 #CH2 bonded to O - oxygen bonded to ring - carbon bonded to O
4 33.0 107.8 #hydrogen from CH3 - carbon from CH3 - hydrogen from CH3
11 50.0 123.2 #benzene carbon - nitrogen 2-xlinks - crosslinked CH2
2 85.0 120.0 #benzene carbon - benzene carbon - carbon from CH3
19 80.0 111.2 #CH ether - crosslinked CH2 - nitrogen 1-xlink
21 35.0 111.0 #benzene carbon - non-xlink nitrogen - hydrogen x-link
30 50.0 118.0 #crosslinked CH2 - nitrogen 2-xlinks - crosslinked CH2
16 85.0 120.0 #benzene carbon - carbon bonded to O - benzene carbon
27 70.0 120.0 #benzene carbon - benzene carbon - nitrogen 1-xlink
28 80.0 109.5 #oxygen before x-link - x-link CH2 - CH ether
13 85.0 120.0 #benzene carbon - benzene carbon - benzene carbon
24 80.0 109.5 #oxygen before x-link - CH ether - x-link CH2
6 70.0 120.0 #graphene carbon - graphene carbon - functional hydroxyl O
33 55.0 108.5 #graphene carbon - functional hydroxyl O - functional hydroxyl H
10 70.0 120.0 #benzene carbon - carbon bonded to O - oxygen bonded to ring
39 70.0 120.0 #graphene carbon - graphene carbon - functional epoxide O
3 70.0 120.0 #benzene carbon - benzene carbon - CH2 bonded to ring
34 55.0 108.5 #CH ether - crosslinked O - crosslinked H
31 80.0 109.5 #CH2 bonded to O - CH ether - oxygen before x-link
9 35.0 120.0 #benzene carbon - benzene carbon - benzene hydrogen
36 50.0 123.2 #benzene carbon - nitrogen 1-xlink - crosslinked CH2
22 70.0 120.0 #benzene carbon - benzene carbon - nitrogen 2-xlinks
35 100.0 111.8 #graphene carbon - functional epoxide O - graphene carbon
18 100.0 111.8 #CH ether - oxygen before x-link - x-link CH2
7 63.0 111.1 #CH2 bonded to O - CH ether - crosslinked CH2

Dihedral Coeffs
dihedraltype K1 K2 K3 K4
31 0.0 0.000 0.000 0.0 #benzene carbon - nitrogen 2-xlinks - crosslinked CH2 - CH ether
40 0.0 2.650 0.000 0.0 #carbon from CH3 - benzene carbon - benzene carbon - nitrogen 1-xlink
25 0.0 0.200 1.450 0.0 #benzene carbon - benzene carbon - carbon bonded to O - oxygen bonded to ring
38 0.0 0.000 0.000 0.0 #CH ether - crosslinked CH2 - nitrogen 2-xlinks - crosslinked CH2
4 0.0 2.650 0.000 0.0 #graphene carbon - graphene carbon - graphene carbon - graphene carbon
29 0.0 1.800 0.000 0.0 #oxygen bonded to ring - CH2 bonded to O - CH ether - crosslinked O
8 0.0 1.000 2.000 0.0 #benzene carbon - carbon bonded to O - oxygen bonded to ring - CH2 bonded to O
7 0.0 2.650 0.000 0.0 #benzene carbon - benzene carbon - benzene carbon - benzene carbon
41 0.0 2.650 0.000 0.0 #benzene carbon - benzene carbon - nitrogen 1-xlink - hydrogen x-link
36 0.0 2.650 0.000 0.0 #benzene carbon - benzene carbon - benzene carbon - nitrogen 2-xlinks
34 0.0 2.650 0.000 0.0 #oxygen bonded to ring - CH2 bonded to O - CH ether - crosslinked CH2
1 0.0 0.200 1.450 0.0 #CH ether - CH2 bonded to O - oxygen bonded to ring - carbon bonded to O
10 0.0 2.650 0.000 0.0 #benzene carbon - benzene carbon - CH2 bonded to ring - benzene carbon
6 0.0 0.200 1.450 0.0 #graphene carbon - graphene carbon - functional epoxide O - graphene carbon
14 0.0 1.800 0.000 0.0 #benzene hydrogen - benzene carbon - carbon bonded to O - oxygen bonded to ring
12 0.0 2.650 0.000 0.0 #benzene hydrogen - benzene carbon - benzene carbon - CH2 bonded to ring
48 0.0 1.800 0.000 0.0 #oxygen before x-link - CH ether - CH2 bonded to O - oxygen bonded to ring
5 0.0 1.000 2.000 0.0 #graphene carbon - graphene carbon - graphene carbon - functional hydroxyl O

79
46 0.0 0.00 1.000 0.0 #CH2 bonded to O - CH ether - crosslinked O - crosslinked H
55 0.0 0.000 0.000 0.0 #benzene carbon - nitrogen 1-xlink - crosslinked CH2 - CH ether
35 0.0 2.650 0.000 0.0 #carbon from CH3 - benzene carbon - benzene carbon - nitrogen 2-xlinks
43 0.0 0.200 1.450 0.0 #CH2 bonded to O - CH ether - oxygen before x-link - x-link CH2
15 0.0 2.650 0.000 0.0 #benzene carbon - benzene carbon - benzene carbon - carbon bonded to O
22 0.0 2.650 0.000 0.0 #benzene hydrogen - benzene carbon - benzene carbon - carbon bonded to O
9 0.0 2.650 0.000 0.0 #oxygen bonded to ring - CH2 bonded to O - CH ether - x-link CH2
2 0.0 1.800 0.000 0.0 #functional hydroxyl O - graphene carbon - graphene carbon - functional hydroxyl O
54 0.0 0.200 1.450 0.0 #graphene carbon - functional epoxide O - graphene carbon - graphene carbon
19 0.0 2.650 0.000 0.0 #benzene carbon - benzene carbon - benzene carbon - benzene hydrogen
44 0.0 0.000 0.924 0.0 #benzene carbon - benzene carbon - carbon from CH3 - hydrogen from CH3
17 0.0 1.000 2.000 0.0 #benzene carbon - carbon bonded to O - benzene carbon - benzene hydrogen
3 0.0 2.650 0.000 0.0 #benzene hydrogen - benzene carbon - benzene carbon - benzene hydrogen
13 0.0 2.650 0.000 0.0 #benzene carbon - carbon bonded to O - benzene carbon - benzene carbon
56 0.0 2.650 0.000 0.0 #CH2 bonded to ring - benzene carbon - benzene carbon - nitrogen 1-xlink
37 0.0 1.000 2.000 0.0 #CH2 bonded to O - CH ether - x-link CH2 - oxygen before x-link
21 0.0 2.650 0.000 0.0 #benzene carbon - benzene carbon - benzene carbon - carbon from CH3
51 0.0 0.00 4.00 0.0 #CH2 bonded to O - CH ether - crosslinked CH2 - nitrogen 1-xlink
18 0.0 0.00 1.000 0.0 #graphene carbon - graphene carbon - functional hydroxyl O - functional hydroxyl H
27 0.0 2.650 0.000 0.0 #benzene carbon - CH2 bonded to ring - benzene carbon - benzene carbon
52 0.0 2.650 0.000 0.0 #benzene carbon - benzene carbon - non-xlink nitrogen - hydrogen x-link
49 0.0 2.650 0.000 0.0 #non-xlink nitrogen - benzene carbon - benzene carbon - carbon from CH3
50 0.0 2.650 0.000 0.0 #benzene carbon - benzene carbon - CH2 bonded to ring - CH3
16 0.0 2.650 0.000 0.0 #benzene carbon - benzene carbon - benzene carbon - CH2 bonded to ring
53 0.0 1.800 0.000 0.0 #functional hydroxyl O - graphene carbon - graphene carbon - functional epoxide O
33 0.0 13.6 0.000 0.0 #benzene carbon - benzene carbon - nitrogen 1-xlink - crosslinked CH2
30 0.0 2.650 0.000 0.0 #benzene carbon - benzene carbon - carbon bonded to O - benzene carbon
32 0.0 2.650 0.000 0.0 #CH2 bonded to ring - benzene carbon - benzene carbon - non-xlink nitrogen
23 0.0 0.00 1.000 0.0 #crosslinked CH2 - CH ether - crosslinked O - crosslinked H
42 0.0 0.00 4.00 0.0 #nitrogen 1-xlink - crosslinked CH2 - CH ether - crosslinked O
20 0.0 0.00 4.00 0.0 #crosslinked O - CH ether - crosslinked CH2 - nitrogen 2-xlinks
47 0.0 2.650 0.000 0.0 #CH2 bonded to ring - benzene carbon - benzene carbon - nitrogen 2-xlinks
11 0.0 1.000 2.000 0.0 #graphene carbon - graphene carbon - graphene carbon - functional epoxide O
39 0.0 2.650 0.000 0.0 #benzene carbon - benzene carbon - benzene carbon - nitrogen 1-xlink
26 0.0 0.00 4.00 0.0 #CH2 bonded to O - CH ether - crosslinked CH2 - nitrogen 2-xlinks
28 0.0 0.000 0.000 0.0 #CH ether - crosslinked CH2 - nitrogen 1-xlink - hydrogen x-link
45 0.0 2.650 0.000 0.0 #benzene carbon - benzene carbon - benzene carbon - non-xlink nitrogen
24 0.0 13.6 0.000 0.0 #benzene carbon - benzene carbon - nitrogen 2-xlinks - crosslinked CH2

Pair Coeffs
epsilon sigma
1 0.0700 3.5500 #benzene carbon
2 0.0300 2.4200 #benzene hydrogen
3 0.180 3.9050 #CH2 bonded to ring
4 0.1180 3.8000 #CH2 bonded to O
5 0.1700 3.8000 #oxygen before x-link
6 0.0800 3.8000 #CH ether
7 0.1700 3.0700 #oxygen bonded to ring
8 0.0700 3.5500 #carbon bonded to O
9 0.1750 3.9050 #CH3
10 0.1700 3.2500 #non-xlink nitrogen
11 0.0000 0.0000 #hydrogen x-link
12 0.0660 3.5000 #carbon from CH3
13 0.0300 2.5000 #hydrogen from CH3
14 0.0700 3.5500 #graphene carbon
80
15 0.1180 3.8000 #x-link CH2
16 0.1700 3.8000 #functional hydroxyl O
17 0.0000 0.0000 #functional hydroxyl H
18 0.1700 3.8000 #functional epoxide O
19 0.1700 3.2500 #nitrogen 1-xlink
20 0.1700 3.8000 #crosslinked O
21 0.1180 3.8000 #crosslinked CH2
22 0.0000 0.0000 #crosslinked H
23 0.1700 3.2500 #nitrogen 2-xlinks

Atoms
id molid type x y z
1 1 3 19.8229454463 96.7510377011 4.0837248078
2 1 1 21.3689651647 96.8215734085 4.26562421574
3 1 1 21.9530165725 98.0891881434 4.4079375841
.....FILE INCLUDES ALL ATOM POSITIONS, BONDS, ANGLES, AND
DIHEDRALS...........

81
 APPENDIX B

Copyright Agreements and Permissions for Figures

B.1. Permission to use Figure 1.2

82
B.2 Copyright Agreement for Text and Figures in Chapter 2

ELSEVIER LICENSE
TERMS AND CONDITIONS

Aug 20, 2014

This is a License Agreement between Cameron M Hadden ("You") and Elsevier

("Elsevier") provided by Copyright Clearance Center ("CCC"). The license consists of
your order details, the terms and conditions provided by Elsevier, and the payment
terms and conditions.
All payments must be made in full to CCC. For payment instructions, please see
information listed at the bottom of this form.
Elsevier Limited
Supplier The Boulevard,Langford Lane
Kidlington,Oxford,OX5 1GB,UK
Registered Company
1982084
Number
Customer name Cameron M Hadden
Customer address 42063 Superior Rd.
HOUGHTON, MI 49931
License number 3452990080935
License date Aug 20, 2014
Licensed content publisher Elsevier
Licensed content
Composites Science and Technology
publication
Molecular modeling of EPON-862/graphite composites:
Licensed content title
Interfacial characteristics for multiple crosslink densities
C.M. Hadden,B.D. Jensen,A. Bandyopadhyay,G.M.
Licensed content author
Odegard,A. Koo,R. Liang
Licensed content date 4 March 2013
Licensed content volume
76
number
Licensed content issue
n/a
number
Number of pages 8
Start Page 92
End Page 99
Type of Use reuse in a thesis/dissertation

83
Portion full article
Format both print and electronic
Are you the author of this
Yes
Elsevier article?
Will you be translating? No
Title of your Molecular Modeling of EPON 862-DETDA / Carbon
thesis/dissertation Composites
Expected completion date Aug 2014
Estimated size (number of
90
pages)
Elsevier VAT number GB 494 6272 12
Permissions price 0.00 USD
VAT/Local Sales Tax 0.00 USD / 0.00 GBP
Total 0.00 USD
Terms and Conditions

84
 INTRODUCTION

1. The publisher for this copyrighted material is Elsevier. By clicking "accept" in

connection with completing this licensing transaction, you agree that the following
terms and conditions apply to this transaction (along with the Billing and Payment
terms and conditions established by Copyright Clearance Center, Inc. ("CCC"), at the
time that you opened your Rightslink account and that are available at any time at
http://myaccount.copyright.com).

GENERAL TERMS

2. Elsevier hereby grants you permission to reproduce the aforementioned material

subject to the terms and conditions indicated.

3. Acknowledgement: If any part of the material to be used (for example, figures) has
appeared in our publication with credit or acknowledgement to another source,
permission must also be sought from that source. If such permission is not obtained
then that material may not be included in your publication/copies. Suitable
acknowledgement to the source must be made, either as a footnote or in a reference list
at the end of your publication, as follows:

“Reprinted from Publication title, Vol /edition number, Author(s), Title of article / title
of chapter, Pages No., Copyright (Year), with permission from Elsevier [OR
APPLICABLE SOCIETY COPYRIGHT OWNER].” Also Lancet special credit -
“Reprinted from The Lancet, Vol. number, Author(s), Title of article, Pages No.,
Copyright (Year), with permission from Elsevier.”

4. Reproduction of this material is confined to the purpose and/or media for which
permission is hereby given.

5. Altering/Modifying Material: Not Permitted. However figures and illustrations may

be altered/adapted minimally to serve your work. Any other abbreviations, additions,
deletions and/or any other alterations shall be made only with prior written
authorization of Elsevier Ltd. (Please contact Elsevier at [email protected])

6. If the permission fee for the requested use of our material is waived in this instance,
please be advised that your future requests for Elsevier materials may attract a fee.

7. Reservation of Rights: Publisher reserves all rights not specifically granted in the
combination of (i) the license details provided by you and accepted in the course of this
licensing transaction, (ii) these terms and conditions and (iii) CCC's Billing and
Payment terms and conditions.

8. License Contingent Upon Payment: While you may exercise the rights licensed
immediately upon issuance of the license at the end of the licensing process for the
85
transaction, provided that you have disclosed complete and accurate details of your
proposed use, no license is finally effective unless and until full payment is received
from you (either by publisher or by CCC) as provided in CCC's Billing and Payment
terms and conditions. If full payment is not received on a timely basis, then any license
preliminarily granted shall be deemed automatically revoked and shall be void as if
never granted. Further, in the event that you breach any of these terms and conditions
or any of CCC's Billing and Payment terms and conditions, the license is automatically
revoked and shall be void as if never granted. Use of materials as described in a
revoked license, as well as any use of the materials beyond the scope of an unrevoked
license, may constitute copyright infringement and publisher reserves the right to take
any and all action to protect its copyright in the materials.

9. Warranties: Publisher makes no representations or warranties with respect to the

licensed material.

10. Indemnity: You hereby indemnify and agree to hold harmless publisher and CCC,
and their respective officers, directors, employees and agents, from and against any and
all claims arising out of your use of the licensed material other than as specifically
authorized pursuant to this license.

11. No Transfer of License: This license is personal to you and may not be sublicensed,
assigned, or transferred by you to any other person without publisher's written
permission.

12. No Amendment Except in Writing: This license may not be amended except in a
writing signed by both parties (or, in the case of publisher, by CCC on publisher's
behalf).

13. Objection to Contrary Terms: Publisher hereby objects to any terms contained in
any purchase order, acknowledgment, check endorsement or other writing prepared by
you, which terms are inconsistent with these terms and conditions or CCC's Billing and
Payment terms and conditions. These terms and conditions, together with CCC's
Billing and Payment terms and conditions (which are incorporated herein), comprise
the entire agreement between you and publisher (and CCC) concerning this licensing
transaction. In the event of any conflict between your obligations established by these
terms and conditions and those established by CCC's Billing and Payment terms and
conditions, these terms and conditions shall control.

14. Revocation: Elsevier or Copyright Clearance Center may deny the permissions
described in this License at their sole discretion, for any reason or no reason, with a full
refund payable to you. Notice of such denial will be made using the contact
information provided by you. Failure to receive such notice will not alter or invalidate
the denial. In no event will Elsevier or Copyright Clearance Center be responsible or
liable for any costs, expenses or damage incurred by you as a result of a denial of your

86
permission request, other than a refund of the amount(s) paid by you to Elsevier and/or
Copyright Clearance Center for denied permissions.

LIMITED LICENSE

The following terms and conditions apply only to specific license types:

15. Translation: This permission is granted for non-exclusive world English rights
only unless your license was granted for translation rights. If you licensed translation
rights you may only translate this content into the languages you requested. A
professional translator must perform all translations and reproduce the content word for
word preserving the integrity of the article. If this license is to re-use 1 or 2 figures then
permission is granted for non-exclusive world rights in all languages.

16. Posting licensed content on any Website: The following terms and conditions
apply as follows: Licensing material from an Elsevier journal: All content posted to the
web site must maintain the copyright information line on the bottom of each image; A
hyper-text must be included to the Homepage of the journal from which you are
licensing at http://www.sciencedirect.com/science/journal/xxxxx or the Elsevier
homepage for books at http://www.elsevier.com; Central Storage: This license does not
include permission for a scanned version of the material to be stored in a central
repository such as that provided by Heron/XanEdu.

Licensing material from an Elsevier book: A hyper-text link must be included to the
Elsevier homepage at http://www.elsevier.com . All content posted to the web site must
maintain the copyright information line on the bottom of each image.

Posting licensed content on Electronic reserve: In addition to the above the

following clauses are applicable: The web site must be password-protected and made
available only to bona fide students registered on a relevant course. This permission is
granted for 1 year only. You may obtain a new license for future website posting.

For journal authors: the following clauses are applicable in addition to the above:
Permission granted is limited to the author accepted manuscript version* of your paper.

*Accepted Author Manuscript (AAM) Definition: An accepted author manuscript

(AAM) is the author’s version of the manuscript of an article that has been accepted for
publication and which may include any author-incorporated changes suggested through
the processes of submission processing, peer review, and editor-author
communications. AAMs do not include other publisher value-added contributions such
as copy-editing, formatting, technical enhancements and (if relevant) pagination.

You are not allowed to download and post the published journal article (whether PDF
or HTML, proof or final version), nor may you scan the printed edition to create an
87
electronic version. A hyper-text must be included to the Homepage of the journal from
which you are licensing at http://www.sciencedirect.com/science/journal/xxxxx. As
part of our normal production process, you will receive an e-mail notice when your
article appears on Elsevier’s online service ScienceDirect (www.sciencedirect.com).
That e-mail will include the article’s Digital Object Identifier (DOI). This number
provides the electronic link to the published article and should be included in the
posting of your personal version. We ask that you wait until you receive this e-mail and
have the DOI to do any posting.

Posting to a repository: Authors may post their AAM immediately to their employer’s
institutional repository for internal use only and may make their manuscript publically
available after the journal-specific embargo period has ended.

Please also refer to Elsevier's Article Posting Policy for further information.

18. For book authors the following clauses are applicable in addition to the
above: Authors are permitted to place a brief summary of their work online only.. You
are not allowed to download and post the published electronic version of your chapter,
nor may you scan the printed edition to create an electronic version. Posting to a
repository: Authors are permitted to post a summary of their chapter only in their
institution’s repository.

20. Thesis/Dissertation: If your license is for use in a thesis/dissertation your thesis

may be submitted to your institution in either print or electronic form. Should your
thesis be published commercially, please reapply for permission. These requirements
include permission for the Library and Archives of Canada to supply single copies, on
demand, of the complete thesis and include permission for UMI to supply single copies,
on demand, of the complete thesis. Should your thesis be published commercially,
please reapply for permission.

Elsevier Open Access Terms and Conditions

Elsevier publishes Open Access articles in both its Open Access journals and via its
Open Access articles option in subscription journals.

Authors publishing in an Open Access journal or who choose to make their article
Open Access in an Elsevier subscription journal select one of the following Creative
Commons user licenses, which define how a reader may reuse their work: Creative
Commons Attribution License (CC BY), Creative Commons Attribution – Non
Commercial - ShareAlike (CC BY NC SA) and Creative Commons Attribution – Non
Commercial – No Derivatives (CC BY NC ND)

Terms & Conditions applicable to all Elsevier Open Access articles:

88
Any reuse of the article must not represent the author as endorsing the adaptation of the
article nor should the article be modified in such a way as to damage the author’s
honour or reputation.

The author(s) must be appropriately credited.

If any part of the material to be used (for example, figures) has appeared in our
publication with credit or acknowledgement to another source it is the responsibility of
the user to ensure their reuse complies with the terms and conditions determined by the
rights holder.

Additional Terms & Conditions applicable to each Creative Commons user

license:

CC BY: You may distribute and copy the article, create extracts, abstracts, and other
revised versions, adaptations or derivative works of or from an article (such as a
translation), to include in a collective work (such as an anthology), to text or data mine
the article, including for commercial purposes without permission from Elsevier

CC BY NC SA: For non-commercial purposes you may distribute and copy the article,
create extracts, abstracts and other revised versions, adaptations or derivative works of
or from an article (such as a translation), to include in a collective work (such as an
anthology), to text and data mine the article and license new adaptations or creations
under identical terms without permission from Elsevier

CC BY NC ND: For non-commercial purposes you may distribute and copy the article
and include it in a collective work (such as an anthology), provided you do not alter or
modify the article, without permission from Elsevier

Any commercial reuse of Open Access articles published with a CC BY NC SA or CC

BY NC ND license requires permission from Elsevier and will be subject to a fee.

Commercial reuse includes:

Promotional purposes (advertising or marketing)

Commercial exploitation ( e.g. a product for sale or loan)

Systematic distribution (for a fee or free of charge)

Please refer to Elsevier's Open Access Policy for further information.

21. Other Conditions:

v1.6
89
You will be invoiced within 48 hours of this transaction date. You may pay your
invoice by credit card upon receipt of the invoice for this transaction. Please
follow instructions provided at that time.
To pay for this transaction now; please remit a copy of this document along with
your payment. Payment should be in the form of a check or money order
referencing your account number and this invoice number RLNK501381820.
Make payments to "COPYRIGHT CLEARANCE CENTER" and send to:

Copyright Clearance Center

Dept 001
P.O. Box 843006
Boston, MA 02284-3006
Please disregard electronic and mailed copies if you remit payment in advance.
Questions? [email protected] or +1-855-239-3415 (toll free in the US)
or +1-978-646-2777.
Gratis licenses (referencing $0 in the Total field) are free. Please retain this
printable license for your reference. No payment is required.

90