A study of thermo-mechanical properties of the cross-linked epoxy: An atomistic simulation

Accepted Manuscript A Study of Thermo-mechanical Properties of The Cross-linked Epoxy: An Atomistic Simulation S. Masoumi, B. Arab, H. Valipour PII:

S0032-3861(15)30057-4

DOI:

10.1016/j.polymer.2015.06.038

Reference:

JPOL 17936

To appear in:

Polymer

Received Date: 13 May 2015 Accepted Date: 16 June 2015

Please cite this article as: Masoumi S, Arab B, Valipour H, A Study of Thermo-mechanical Properties of The Cross-linked Epoxy: An Atomistic Simulation, Polymer (2015), doi: 10.1016/j.polymer.2015.06.038. This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT

A Study of Thermo-mechanical Properties of The Cross-linked Epoxy: An Atomistic Simulation S. Masoumi1,a, B. Arabb, H. Valipoura a

Centre for Infrastructure Engineering and Safety (CIES), School of Civil and Environmental Engineering, UNSW Australia, UNSW Sydney, NSW 2052, Australia.

b

RI PT

Department of Mechanical Engineering, Faculty of Engineering, North Tehran Branch, Islamic Azad University, Tehran, Iran.

SC

Abstract

TE D

M AN U

In this paper, crosslinking process and thermo-mechanical properties of epoxy materials containing diglycidyl ether bisphenol-A (DGEBA) as the epoxy resin and JEFFAMINE® D-230 polyoxypropylenediamine as the hardener have been studied using molecular dynamics simulation. An algorithm to create the crosslinked epoxy has been developed, and the crosslinking process is monitored. The simulation results disclose the effectiveness of the using annealing in the crosslinking process to achieve higher crosslinking densities within shorter reaction radiuses while producing more desirable conformation. The evolution of the complex network of cured epoxy is studied by radial distribution function (RDF). The important properties of the epoxy material such as glass transition temperature (Tg), the coefficient of thermal expansion (CTE), and elastic constants are calculated. Furthermore, the variation of the properties through the evolution process is considered, and the improvement of them is captured. The results obtained from our simulation are in good agreement with the experimental results, revealing the strength of molecular dynamics to predict the material properties. Keywords: Atomistic simulation, Crosslinked Epoxy, Thermo-mechanical properties

EP

1. Introduction

AC C

Epoxy is the product of the reaction between the epoxy resin and a curing agent. Epoxy resins are the class of polymers containing epoxide groups, and they may react with each other or with other curing agents to create an epoxy material with outstanding properties. The curing agents usually include co-reactant groups such as amines and phenols. The result of the reaction between epoxide groups and co-reactant groups in curing agents usually is a thermosetting polymer with superior thermos-mechanical properties. Due to their excellent properties, epoxies have a wide range of applications in various industries including but not limited to aircraft structures, FRP-concrete bonding, electrical components, and shipbuilding. Since the introduction of molecular dynamics (MD) by Alder and Wainwright [1] in the late 50’s and later by Rahman[2], there are numerous research works focusing on the application and improvement of MD, mostly in the field of science. MD is computationally demanding, and it was not until recently that modern computers made it possible for researchers to use its privileges in the 1

Corresponding author 1

ACCEPTED MANUSCRIPT engineering problems such as structural and material analysis. Currently, the crosslinking process of polymers is one of the exciting research areas where MD can be applied to obtain the results that would help the researchers to understand underlying reasons for specific behaviors observed in the experiments. The available results show a remarkable potential for MD method to be considered as a reliable and promising method to predict thermo-mechanical properties of the polymers.

SC

RI PT

The crosslinking process using MD simulation usually follows the process happens naturally in the experiments, where the creation of bonds in shorter reaction radiuses occurs initially. Yarovsky and Evans [3] proposed an algorithm to simulate the crosslinking process using molecular dynamics approach. They have used COMPASS force field to define the interactions between atoms, and a cutoff radius of 6 Å is suggested to create the bonds. They have shown how the volume shrinks with respect to the curing ratio. Also, they have investigated the barrier properties of the epoxy against the penetrant molecules and shown how diffusivity is affected by the interaction between the penetrant molecule and the epoxy structure. Adhesion properties of the epoxy are also examined by means of free energy concept.

AC C

EP

TE D

M AN U

Wu and Xu [4] proposed a crosslinking process to create bonds between epoxy resin and curing agent using molecular dynamics approach. Their results include mechanical properties of the resultant epoxy material. In another work, Wu and Xu [5] has investigated the epoxy structure using atomistic simulation. They have used COMPASS force field to work out the properties of the epoxy. To overcome the difficulty of obtaining Tg with their model, they have used experimental data to initiate their simulation in an NVT-MD ensemble. They studied the structure of the material using radial distribution function (RDF) as well as mean squared displacement (MSD) and vector auto correction function. Tack and Ford [6] considered the epoxy as a mixture of DGEBF and DETA. They calculated the mechanical properties of the epoxy using MD simulation. They used two different forcefields COMPASS and cff91 to examine the efficiency of each one. It is shown that COMPASS can predict the properties more closely to experimental results than cff91. Varshney et. al. [7] proposed four approaches to create the bonds between epoxy resin EPON-826 and curing agent DETDA based on the reactivity of the primary and secondary amine groups. They investigated the thermal and mechanical properties of the system using CVFF forcefield. Their results show a good agreement with experimental results. Nouri and Ziaei-Rad [8] studied the properties of crosslinked epoxy using a simple force-field expression. They have shown their method is fast and efficient to calculate thermo-mechanical properties. They have investigated the structure of their network using fractal dimension. Dynamics method can be used to calculate the mechanical properties of the materials. Dynamics method is capable of capturing the temperature dependency of the properties. Yang and Qu [9] studied thermo-mechanical properties of crosslinked epoxy by atomistic simulation. They have used PCFF forcefield to define the interactions between atoms. RDF has been used to study the structure of the system. They have also investigated the dependence of the thermos-mechanical properties to the crosslinking density as well as temperature. Sirk et al. [10] have investigated molecular dynamics simulation as well as experimental observation of the high strain rate mechanical properties of a cross-linked epoxy by mixing 972 molecules of epoxy resin and 486 molecules of hardener. They used AMBER force-field to define the interactions between atoms. Glass transition temperature obtained from their simulation was higher than the experimental results. They have used WLF equation [11] to relate the difference between the Tg from simulation and experiments to the 2

ACCEPTED MANUSCRIPT

M AN U

SC

RI PT

cooling rate. The shift factor from WLF equation predicted that Tg in the simulation must be higher than the experimental results. They compared the Tg modified with the shift factor to experimental results, and good agreement has been reported using this approach. The WLF equation is based on the experimental observation and is suggested to be useful to compare the results from the simulation and experiments[12]. They have also worked out mechanical properties of the epoxy with various strain rates. Shokuhfar and Arab [13] used a 4 step crosslinking procedure to build the epoxy network with different crosslinking densities. They have calculated the mechanical properties of the systems using the static method. The dependency of the properties to the crosslinking density is considered. Their results were close to those obtained from experiments. They have studied the structure of the epoxy system using RDF. Monk et. al. [14] have established a sensitivity analysis of the crosslinking phenolic resins using molecular dynamics approach. They have considered different effective parameters including time, initial volume, crosslinking approach, equilibration temperature, and the ensemble while using full factorial analysis to measure the effect of each parameter. Their results suggest the importance of each parameter for the calculation of the material properties. Soni et al. [15] studied the effect of the chain length of hardeners on the thermal properties of the epoxy materials. Their molecular structure included 150 epoxy resins and 75 hardeners. They found out that the density of the molecular structure would be lower for the longer hardener’s chain while their coefficient of thermal expansion (CTE), as well as glass transition temperature (Tg) decreases. There are more researchs focusing on the properties of the crosslinked epoxy using molecular dynamics [16]–[23]. Reinforcing the crosslinking epoxy with nano-particles has also been considered in other works [24]–[28].

AC C

EP

TE D

In this paper, an efficient algorithm has been developed to create the crosslinked epoxy. The focus of the algorithm is to achieve the highest crosslinking ratio in the shorter reaction radius. This algorithm can help to copy the process usually occurs in reality. The details of the crosslinking are discussed, and the role of annealing process is divulged. The local structure of the new formed network is examined using radial distribution function (RDF). The important properties of the epoxy materials such as glass transition temperature (Tg), the coefficient of thermal expansion (CTE), and young modulus are calculated using the molecular dynamics method. Also, the dependency of the properties to the crosslinking density is investigated. Despite using small size molecular dynamics simulation, by focusing on the reliable crosslinking algorithm and precise simulation techniques available at present, good results are obtained. The outcomes of the simultion are verified with the experimental observations.

2. Crosslinking methodology The crosslinking process in the epoxy material is the process that usually the reaction occurs between the end carbons in the epoxy resin containing epoxide groups and a hardener often containing active amine groups at their both ends. In our simulations, the epoxy resin is diglycidyl ether bisphenol-A (DGEBA) containing two epoxide groups at two ends (Fig.1 a) and Jeffamine® D230 (Fig. 2b), a trademark of Huntsman Corporation. D-230 molecules have repeated oxypropylene units as their backbones, and the amine group has been located on the secondary carbon atoms at each end of the aliphatic polyether chain. In our simulation the repeated number of oxypropylene is considered three. The reaction happens because the end carbon in the epoxy ring has slightly 3

ACCEPTED MANUSCRIPT

M AN U

a

SC

RI PT

positive charge while nitrogen in the amine group has slightly negative charge; the nitrogen would absorb the end carbon while losing a hydrogen atom. The hydrogen atom that has slightly positive charge would be absorbed by the oxygen in the opened epoxy ring. The creation of OH group in the system would contribute to some of the superior properties of the epoxy such as adhesivity. The nitrogen in amine group can react with two carbons at once resulting in the connection between three molecules. Usually, the end carbon reacts with the nitrogen as it is more accessible (Fig. 2). This linking process would lead to a creation of the system that consists of a complex network in which almost all molecules are included. The superior properties of the epoxy materials are because of the existence of such network.

TE D

b

AC C

EP

Fig. 1 Schematic representation of a) DGEBA b) D-230 molecules (x=2.6).

Fig. 2 The ring opening and bond creation process between DGEBA and D-230 molecules.

The interactions between atoms have been determined using well-defined force-field COMPASS, which is an ab initio force-field [29]–[32]. This forcefield stresses more on the accuracy and, therefore, is suitable to compare the simulation results with experimental outcomes. The potential function has two categories: valance terms including diagonal and off-diagonal cross-coupling terms

4

ACCEPTED MANUSCRIPT and non-bond terms including van der Waals (vdW) term described by Leonard-Jones (LJ) 9-6 and a Columbic function to define the electrostatic interactions. The force-field function is given as follow:

[

Etotal = ∑ K 2 (b − bo ) 2 + K 3 (b − bo ) 3 + K 4 (b − bo ) 4

[

+ ∑ H 2 (θ − θ o ) 2 + H 3 (θ − θ o ) 3 + H 4 (θ − θ o ) 4 θ

]

]

+ ∑ [V1 (1 − cos φ ) + V2 (1 − cos 2φ ) + V3 (1 − cos 3φ )] φ

RI PT

b

+ ∑ K x ( x − xo ) 2 + ∑ K bb ' (b − bo )(b'−b'o ) + ∑ Kθθ ' (θ − θ o )(θ '−θ 'o ) x

θ ,θ '

b ,b

+ ∑ K bθ (b − bo )(θ − θ o ) + ∑ (b − bo )[G1 cos φ + G2 cos 2φ + G3 cos 3φ ] b ,θ

b ,φ

SC

+ ∑ (θ − θ o )[F1 cos φ + F2 cos 2φ + F3 cos 3φ ] θ ,φ

∑ K θ (θ − θ

θ ,θ ,φ

+∑ i, j

b

qi q j rij

o

)(θ '−θ 'o )(φ − φo )

 rijo rijo  + ∑ ε ij 2( ) 9 − 3( ) 6  rij  i, j  rij

M AN U

+

(1)

TE D

Where the first fourth terms represents the energy associated with bond (b), angle (ϴ), torsion (φ), and Wilson out-of-plane (x) internal coordinates. The next fifth terms are the cross-coupled terms and the last two terms express the non-bond interaction between atoms. The second to last term is relatively long range Columbic potential defining electrostatic interaction, while the last term is Leonard Jones (LJ 9-6) potential to defining van der Waals interaction.

EP

In this report, Ewald summation method has been used to sum the long-range electrostatic interaction terms while the direct atom-based method is chosen for the short-range vdW interactions. The periodic boundary condition is considered to define the interaction of atoms in boundaries with the neighbor cells.

AC C

2.1. Crosslinking algorithm

All simulations are carried out using Materials Studio software[33]. In order to design an atomistic approach to model the crosslinking process that usually occurs in the reality the following steps has been taken to prepare the structure for the process: Step one: DGEBA and D-230 molecules minimized and then they have been randomly spread into a cell with the ratio of 2:1, respectively. This ratio is because every nitrogen atom is enabling to react with two ending carbons in the resin molecule. At this stage, 100 different configurations has been created and minimized, and the configuration with minimum total energy is chosen as the acceptable conformation. This step is necessary to make sure the molecular structure’s energy is close to the global minimum energy. Amorphous Cell module has been used to mix the molecules into the cell.

5

ACCEPTED MANUSCRIPT Step two: The molecular structure undergoes a short NVT-MD for 50 ps and then an equilibration using NPT-MD for 200 ps in 300K. Then, the annealing process applied to the system to overcome the energy barriers and obtain a total energy that is closer to the structure with global minimum energy. The cooling rate at this stage is 10K/10ps, and it heats up to 550K and cools down to 10K during four cycles. The time step in preparation and the crosslinking process is 1 fs.

AC C

EP

TE D

M AN U

SC

RI PT

The crosslinking process would start when the initial structure is ready. In this report, DGEBA with closed ring has been mixed with hardeners. Whenever the reactants from amine group and epoxide ring come close to each other, epoxy ring would be opened, and the bond would be created to imitate the process that happens in reality (Fig. 2). The flowchart in Fig. 3 shows the algorithm for construction and equilibration of the crosslinking epoxies.

Fig. 3 The crosslinking algorithm.

6

ACCEPTED MANUSCRIPT

TE D

M AN U

SC

RI PT

As it can be seen in the flowchart, first the initial values are set. The minimum and maximum cut-off distance for creating the bond between carbon in epoxide group and nitrogen in the amine group is 4 Å and 13 Å, respectively. The number of iterations in each radial distance is 5. The cut-off distance for short-range vdW type interaction terms in MD simulation is 12.5 Å while Ewald summation method with the accuracy of 0.001 and buffer width of 0.5 Å is considered for the long-range electrostatic interaction. In order to control the temperature, Nose –Hoover-Langevin [34], [35] thermostat is used while Berendsen barostat [36] is applied to control the pressure for NPT-MD simulation. In addition, the reaction between nitrogen with two carbon atoms in epoxide groups has been considered in this algorithm as well. Once the initial parameters are set, the crosslinking process starts. The process in the simulation follows the main parts of the process occurs in the reality. Whenever the distance between nitrogen in amine group and the end carbon in the epoxy ring is less than the cut-off value specified in the algorithm, the epoxy ring would be opened and the bond would be created between the nitrogen and the carbon. At the same time, oxygen would absorb the hydrogen from the amine group. Because the new bond would result in the new formation of the structure, equilibrating is necessary. First the geometry of the structure is optimized to bring the bond length and the structure within a reasonable range; then, a short NVTMD for 20 ps is applied on the system. The new structure might experience energy barriers to taking an optimized shape. Therefore, annealing has been undertaken to overcome this issue. The annealing process at this stage has the cooling rate of 50K/50ps between 550K and 200K. After annealing the structure equilibrated using NPT-MD for 60 ps, making it ready for creating new bonds. This process iterates until either the maximum number of iteration is reached, or the algorithm cannot find any possible pairs for bond creation; then, the cut-off distance increases with a step size of 1 Å. The algorithm would stop the process when the maximum cut-off length is reached or when the desired crosslinking ratio is obtained.

AC C

EP

Fig. 4 shows how the hardener’s molecules types based on the number of the reactions with the carbon atoms emerge and vanish. While the analysis in this report is based on the system including 24 DGEBA molecules and 12 D-230 molecules, for the sake of the more smooth graphs of evolution process of the molecular structure, a system of 48 DGEBA molecules and 24 D-230 molecules were crosslinked using the algorithm explained in the previous section. As it can be seen from the graph, at the beginning of the crosslinking process, there are only 25% of the hardeners in the system without any reaction. The percentage of the hardeners with only one and two reaction is around 55% and 20% contributing to 13.75% and 10% of the total crosslinking density, respectively. The percentage of the hardeners with no reaction falls gradually as the reaction radius and iterations number increases until all crosslinker molecules have at least one reaction at the radius of 6 Å. Interestingly, the percent of the crosslinkers with four reactions starts from the first reaction radius after five iteration, while there are still crosslinkers without any reactions in the system. The percent of the hardener molecules with two reactions shows no particular pattern as some of the molecules with one reaction would join this group while others leave the group to join to the molecules with three reactions. Molecules with three reactions show an increase up to the reaction radius of 7 Å and then decrease consistently with the increase in the reaction radius and the iterations. The decrease in the number of crosslinkers with three reactions explains the increase in the percentage of the crosslinkers with four reactions. As it was expected, the increase in the reaction radius has a profound effect on the variation of the crosslinking percent. However, it is also important to note 7

ACCEPTED MANUSCRIPT that after the creation of the bonds and equilibration, there are still potential reactions in the system. These can significantly change the final molecular conformation and, therefore, cannot be ignored.

100

70 60 50 40

SC

Percent of involvement

80

RI PT

Crosslinkers with 0 reactions Crosslinkers with 1 reactions Crosslinkers with 2 reactions Crosslinkers with 3 reactions Crosslinkers with 4 reactions

90

30

10 0 0

4-1

4-2

4-3

4-4

4-5

5-1

5-2

M AN U

20

5-3

5-4

6-1

6-2

7-1

7-2

7-3

7-4

8-1

9-1

9-2 10-1 13-1

Reaction radious (Å)-Iteration number

Fig. 4 Contribution of the different crosslinkers types to the total conversion of the epoxy at various stages of evolution.

AC C

EP

TE D

Fig. 5 illustrates the evolution of the crosslinked epoxy with and without considering the annealing process in the simulation. It is interesting to note that the shifts in the system are more pronounced at the new reaction radiuses. While the major shifts happen in the new radiuses, there are minor shifts taking place in the same reaction radius but different iteration number. The small shifts are much bigger in the simulation with annealing compare to those in the simulation without annealing, indicating a more smooth change of the crosslinking density and better formation of the molecular structure. In addition, as it can be seen from the graph, at the reaction radius of 9 Å and 10 Å almost 93% and 96% of the reactions have taken place in the system with annealing, respectively. These conversions do not occur until at the reaction radius of 11 Å and 13 Å for the algorithm without annealing. After the reaction radius of 10 Å, the system with annealing has difficulty to find the new potential reactants as 96% of the free sites are already occupied in the system. The final crosslinking density of 98% occurs at the reaction radius of 13 Å. Whereas, the system without annealing finds new potential reactants as the reaction radius increases until it ends up with the final crosslinking density of 96% at the reaction radius of 13 Å. Although the final conversion density is approximately same for the systems with and without annealing by considering the long reaction radius, there is a remarkable distinction between the paths they converge to the final structure. For the annealing case, the reactions take place at shorter reaction radiuses, and more reactions introduced into the system within the same reaction radius in each iteration. It seems the annealing process can imitate the process happens in reality better. Referring to the Fig. 5, and considering the computational cost of crosslinking process and the high strain produced in the system for long reaction radiuses, the choice of final reaction radius of 8-10 Å for the system with annealing seems reasonable and gives a conversion density that naturally happens in the experiments[19]. It should be mentioned that we 8

ACCEPTED MANUSCRIPT could obtain a crosslinking of 100% for a system with 24 DGEBA and 12 D-230 molecules using the algorithm with annealing with the maximum reaction redius of 13 Å. This system is used for all of the analysis in later sections. 100

Cured without annealing Cured with annealing

RI PT

80

70

60

50

40

SC

Percent of total crosslinking

90

30

4-1 4-2 4-3 4-4 4-5 5-1 5-2 5-3 5-4 6-1 6-2 7-1 7-2 7-3 7-4 8-1 9-1 9-2 10-1 11-1 11-2 12-1 13-1 13-2

M AN U

20 0

Reaction radious (Å)-Iteration number

Fig. 5 Evolution of the crosslinking process with and without considering annealing.

3. Local structure: RDF

N

TE D

Pairwise correlation function (Radial Distribution Function (RDT)) is used to study the local structure of the epoxy and its evolution during the curing process. Given atom A is located in the center of a spherical shell with infinitesimal thickness, RDF measures the probability of finding atom B with its center located inside the shell thickness in the distance r. RDF can be formulated as follow:

1 Nα β xα xβ ρ g (r ) = ∑∑ 'δ (r − r j + ri ) N i =1 j =1

AC C

EP

(2)

Where xi denotes the molar fraction of atoms with chemical type i, Ni is total number of atoms of chemical type i, and prime indicates the exclusion of the terms where i=j. RDF would converge to unity as the radius of the sphere increases while it is zero within the radii of the atom at the origin. There are 6 different crosslinking densities of the epoxy, ranging from 0% to 100%. In this report, CR refers to the carbon in the epoxide group in resin which would react with nitrogen from amine group in the curing agent. Fig.6 shows the interaction between O and H in the system. There are two different O play roles in the graph; one is the O in the epoxide group, and the other is O in the curing agent chain. As it can be seen, there are two major peaks in the graph. The first peak around 0.96 Å (O-H bond length) corresponds to the direct bond between O in epoxide group and H after the ring is opened due to the reaction between C and N. By increasing the crosslinking percentage, the probability of finding H 9

ACCEPTED MANUSCRIPT directly bonded to O raises since there are more reactions in the system. The second peak shows the interactions between O and the H connected to adjacent C in the O-C-H sequence. As the crosslinking increase, the second peak falls since the O in the epoxide group lose its contribution to the second peak.

RI PT

Fig.7 illustrates the interactions of N in amine group and H. The first peak around 1.01 Å shows the direct bond between N and H before a reaction takes place. The more N reacts with C, the lower the probability of finding N-H bond would be. The second peak shows the interaction between N and H in the N-C-H sequence in the curing agent and N-CR-H sequence in the reacted molecules.

SC

The epoxide ring-opening process can be captured through the investigation of the interaction between CR and O. As it has been depicted in the Fig. 8, there are two peaks in the graphs; the first peaks around 1.43 Å shows the direct bond between CR and O in the epoxide ring before opening, while the second peak indicates this correlation after the ring opening.

TE D

M AN U

To capture the evolution of the system and the effects of the N on it, the correlation between CR and CR has been depicted in the Fig.9. As it can be seen, for the crosslinking percentages less than 20 percent the probability of finding two CR close together is zero. It means that there is no N in the system reacting with two CR. In 40% crosslinking the probability for both epoxies is nonzero, indicating the presence of the N in the system with two reactions. As this ratio is less than 50%, there are some Ns in the system without any reaction. The difference in the peaks value for ratios less than 87% is because of the crosslinking process in the systems, and no conclusion can be made. However, for the 100% crosslinking system where all nitrogens in the system are most likely to have reacted with two carbon atoms. 18 16 14

EP

15

10

10

8

5

AC C

g(r)

12

0% Crosslinking 20% 40% 56% 87% 100%

6

0 0.9

0.95

1

1.05

4 2 0

0.6

0.7

0.8

0.9

1

r(Å) Fig. 6 Pair functions between O-H for various crosslinking densities.

10

1.1

1.2

1.3

1.4

ACCEPTED MANUSCRIPT

40

40

35

30

30

20

0% Crosslinking 20% 40% 56% 87% 100%

10

25

0 0.9

20

0.95

1

1.05

1.1

15 10 5 0

0

0.5

1

1.5

2

2.5

3

60

60 40

50

20 0 1.3

g(r)

40

20

0

EP

10

0

1.35

1.4

1.45

0.5

1

1.5

2

4.5

5

1.5

1.55

0% Crosslinking 20% 40% 56% 87% 100%

1.6

10

TE D

30

4

M AN U

Fig. 7 Pair functions between N-H for various crosslinking densities.

70

3.5

SC

r(Å)

RI PT

g(r)

45

5

0 2.2

2.5

r(Å)

AC C

Fig. 8 Pair functions between CR-O for various crosslinking densities.

11

2.4

3

2.6

3.5

4

4.5

5

ACCEPTED MANUSCRIPT 30 0% Crosslinking 20% 40% 56% 87% 100%

25

15

10

5

0

0.5

1

1.5

2

2.5

r(Å)

3.5

4

4.5

5

M AN U

Fig. 9 Pair functions between CR-CR for various crosslinking densities.

3

SC

0

RI PT

g(r)

20

4. The glass transition temperature and the coefficient of thermal expansion

TE D

Glass transition temperature is a second-order thermodynamic transition, which can be defined as a temperature range that a thermosetting polymer structure changes from a hard form (glass state) to a more flexible and mobile form (rubber state). Despite the fact that Tg is a temperature range, it is common to report one point, referencing to Tg. This point is usually within the temperature range where the slopes of the curves in the glassy state and rubbery state in the volume-temperature graph meet each other.

AC C

EP

The output of the algorithm from the section 2 is equilibrated using NPT-MD with a modified version of Nose-Hoover thermostat, Nose-Hoover-Langevin themostat [34], [35] with Q ratio of 0.01. The pressure is controlled using Berendsen barostat. The process of cooling system starts after equilibration with different barostat. The barostat developed by Martyna et al [37] has been carried out to control the pressure of the molecular structure in the current simulation. The algorithm allows the length and angels (volume) of the cell vary in the response to external pressure to create the desired pressure. This notion follows the approach outlined in the barostat developed by Parrinello and Rahman [38]. This barostat introduces more degree of freedom to the system and would lead to more natural simulation especially in the case of solid materials simulation where the cell shape matters. We found this barostat more suitable for our cooling process. In this report, the temperature interval of 10K has been chosen while the system equilibrated for 300 ps NPT-MD. The trend of the average density in the system is the criteria here to whether accept the results or reequilibrating it. If the average density of the system remains constant, the density would be recorded; otherwise, the system would go just through another 300ps NPT-MD. Fig. 10 shows the variation of the specific volume versus the temperature for four different crosslinking densities of DGEBA/D-230. As mentioned before, Tg can be determined using specific volume versus temperature diagram by finding the point where the slope of the curves in rubbery

12

ACCEPTED MANUSCRIPT and glass states meet each other. Two lines with the least standard deviation from the average have been used to fit the data points for each crosslinking ratio, and the intersection of them would identify Tg.

1.1

1

0% Crosslinked 40% 87% 100% Fitted line

TE D

3

0.95

0.9

0.85

300

350

400

450

500

550

Temperature (K)

AC C

250

EP

Specific volume (cm /g)

1.05

M AN U

SC

RI PT

The increase in the crosslinking density leads to a higher Tg (Table. 1). Since the transition is believed to happen as a result of a large-scale movement of the polymer segments at a critical fractional free volume [39] and also considering the fact that the more crosslinking density means the larger polymer segment, it is not surprising to observe the higher Tg for higher crosslinking densities. The Tg predicted for DGEBA/D-230 epoxy with 100% and 87% crosslinking density is around 383K and 373K, respectively. Considering the fact that the crosslinking density of 100% is less likely to happen in the experiments due to the topological constraints, the results from the simulation for 87% crosslinking density is outstanding, while the conformation of 100% crosslinking is in good agreement with the experimental results as outlined in the Table. 1. The higher value anticipated by the simulation compared to the experimental results is reported by other researchers as well. The higher value of the Tg predicted by MD can be attributed to the perfect modeling of the molecular structure, the mechanism of the temperature change in the simulation, and cooling rate dependency of the results. More research is still required in this area.

Fig. 10 specific volume versus temperature for various crosslinking densities.

The coefficient of volumetric thermal expansion (CTE) is the amount of change in size of a matter in the response to temperature variation, and it can be described by the following equation:

αv =

1  ∂v    v0  ∂T  p

(3)

13

ACCEPTED MANUSCRIPT

Table 1 Glass transition temperature and CTE for different crosslinking densities.

Material description

Tg (K)

0% Crosslinked 40% Crosslinked 87% Crosslinked 100% Crosslinked Experiments [40] [41] [42] [43] MD simulations 2 [10] 3 [15]

345 361 373 384

αv Glassy

Rubbery 7.64 7.34 4.93 4.53

M AN U

SC

3.24 3.82 1.92 1.83

RI PT

where v0 is the specific volume at 300K in this report. The subscript p denotes the constant pressure in the system while differentiating the volume respect to the temperature. In Table 1, the CTE for the glassy and rubbery states is provided. The results show a general trend of decrease in the CTE for the glassy and rubbery states as the crosslinking density increases. As it can be seen the CTE in the rubbery state is more than two times higher than CTE in glassy indicating soften molecular structure in the rubbery state. The trend of the results is similar to those obtained by Bandyopadhyay [16]. The absolute value of the CTE is in good agreement with experimental results.

360 353 1 338-363 360

1.65 -

6.15 -

445 385

1.95 1.5

5 3.66

1

The range is based on the different amine curative required to react with 100 pbw of epoxy resin (phr). The simulation cell included 972 DGEBA molecules and 486 D-230 molecules with the crosslinking density of 100%. Amber force-field has been used for the analysis. 3 The simulation cell included 150 DGEBA molecules and 75 D-230 molecules with the crosslinking density of 100%. Amber force-field has been used for the analysis.

TE D

2

EP

5. Mechanical properties

AC C

In molecular dynamics, there are mainly three approaches to work out mechanical properties of the materials namely statics, dynamics, and fluctuations method. In the static approach, the contribution of the change of entropy as well as vibrational movements of atoms are assumed to be negligible on the elastic constant calculation. The results obtained from this method are in good agreements with experimental results [44]. On the other hand, dynamics approach takes into account the contribution of vibrational frequencies ad entropy change in the calculation of mechanical properties. Therefore, it can be used to calculate the mechanical properties at various temperatures. However, dynamics method needs long-duration simulations [45], and there are strain fluctuations which result in the uncertainty of the outputs. The third way is the fluctuation method, which can be used to calculate the properties obtained in dynamics method. This process is computationally expensive though it can be improved using modified fluctuation approaches. In this paper, the static method has been carried out to work out the mechanical properties of the epoxy material.

14

ACCEPTED MANUSCRIPT

 N   T  T  ∑ mi ( v i v i )  +  ∑ rij f ij    i< j   i =1

(4)

SC

1 σ=− V0

RI PT

The molecular structure configuration coming from the protocol explained in section 1 is considered to compute the mechanical properties. The system is first undergone an NPT-MD equilibration until the thermodynamics variable states fluctuate about a constant value. The data is recorded in every 5ps and the properties are calculated over the last 100ps of the simulation. For equilibration MD, the Ewald summation with the accuracy of 10-4 is considered to sum the long-ranged electrostatic terms. The cut-off distance of 15.5 Å is considered to sum the vdW terms. The Nose-Hoover-Langevin method is chosen to control the temperature, and the Berendsen barostat is applied to adjust the pressure in the system. The stresses in the system can be obtained by so-called virial method [46] as follow:

C mnpq =

∂σ mn ∂ε pq

= T ,e pq

1 ∂2 A V0 ∂ε mn ∂ε pq

T ,ε mn ,ε pq

M AN U

where N is the number of atoms in the system; mi and vi are the mass and velocity of atom I; fij and rij denote the force and distance between atom i and j. V0 is the undeformed system volume. Once the stresses in the system are determined, the elastic stiffness coefficients can be worked out as follow:

(5)

Cij =

TE D

where A is the Helmholtz free energy. Taking into account the symmetry of the stress and strain tensors, for isotropic materials the stiffness tensor becomes a 6×6 matrix. Applying the static method assumptions, the stiffness matrix can be expressed as:

1 ∂ 2U V0 ∂ε i ∂ε j

(6)

λ λ 0 λ + 2µ λ 0 λ λ + 2µ 0 µ 0 0

AC C

λ + 2 µ λ  λ  0 0  0

EP

Assuming the molecular structure is behaving like an isotropic material, the stiffness matrix can be defined using Lamé constants as follow:

0 0

0 0

0 0

0 0 0  0 0  0 0 µ 0  0 µ  0

(7)

In the equation 7, λ and μ are Lamé constants and are defined as:

1 6 1 µ = (C11 + C 22 + C33 ) 3

λ = (C12 + C13 + C21 + C23 + C31 + C32 )

(8a) (8b) 15

ACCEPTED MANUSCRIPT Finally, the elastic modulus, Poisson’s ratio and shear modulus can be computed using Lamé constants as follow:

µ (3λ + 2µ ) λ+µ λ υ= 2(λ + µ ) G=µ E=

(9a)

RI PT

(9b) (9c)

M AN U

SC

Assuming the material behave as an isotropic material is necessary for our simulation as the size of the system is small, and the stiffness matrix might not show exactly same material properties in each direction. Therefore, the Lamé constants are used to calculate the average properties of different directions. constant strain method is used to obtain the mechanical properties. The molecular structure is subjected to twelve deformations including three pairs in the uniaxial direction (compression and tension) and three pairs involving pure shear.

TE D

The density and the elastic constants are given in Table 2. The reports from the experiments and other MD simulations are all at 300 K temperature. As it can be seen, the density of the molecular structure raises as the crosslinking density increases. The mechanical properties show consistent improvement during the crosslinking process. The static method used in here ignores the effect of temperature on the total energy of the system, however, the structure has been equilibrated at this temperature and therefore it is reasonable to compare the results in this temperature. The absolute values of the mechanical properties are in a good agreement with those obtained from the experiments. The overestimation of the Young’s Modulus in MD simulations is reported by other researchers as well and can be attributed to the perfect modeling of molecular structure that is not the exact modeling of the real materials. Furthermore, ignoring the dynamics effect on the total energy of the system while calculating the properties.

EP

Table 2 The mechanical properties of various crosslinking densities.

Density (g/cm3)

0% crosslinked 40% crosslinked 87% crosslinked 100% crosslinked Experiments [43] [41] [42] DMA [10] MD Simulations [10]3 [15]

AC C

Material description

Poisson ratio

1.103 1.115 1.129 1.138

Elastic modulus (GPa) 2.56 3.09 3.7 3.9

1.29 -

2.641 3.1±0.2 2.88-3.572 2.8

-

1.121 1.124

3.9-5.5 -

0.33-0.41 -

1

The value is extracted from the linear extrapolation of the results. The Same trend as mentioned in Table 1 caption. 3 The values obtained from dynamics method with different strain rates. 2

16

0.32 0.32 0.35 0.3

ACCEPTED MANUSCRIPT 6. Conclusion

References

TE D

M AN U

SC

RI PT

This work reports the molecular dynamics simulation of the crosslinked epoxy material. The epoxy in this report is the mixture of 24 DGEBA molecules and 12 D-230 hardener molecules. This report shows the ability of molecular dynamics simulation as a reliable method for predicting the material properties. A crosslinking algorithm is developed, and the details are explained. The effect of annealing in the system is examined and it is found out that the annealing process would result in more reactions in the same radius in each iteration; also, the total conversion is becomes higher in shorter reaction radiuses compare to the algorithm without annealing. The emergence and disappearance of the different crosslinker types also investigated. The evolution of the local structure of the crosslinked epoxy is studied using radial distribution function (RDF). Glass transition temperature and coefficient of volumetric thermal expansion were calculated for different crosslinking densities using a barostat which let the shape of cell to be changed to control the pressure. The approach showed an excellent convergence and the results were in good agreement with experimental observations. The static method is used to work out the elastic constants, and the results were close to the experimental results. In all cases, the improvement of the thermomechanical properties is observed as the curing process continuous. The focus of this study was the use of more precise simulation techniques for the small size system. The choice of force-field and the equilibration algorithms had a profound influence on the outstanding results obtained from molecular dynamics simulation of the small size system in this report.

B. J. Alder and T. E. Wainwright, “Phase Transition for a Hard Sphere System,” J. Chem. Phys., vol. 27, no. 5, p. 1208, 1957.

[2]

A. Rahman, “Correlations in the motion of atoms in liquid silicon,” Phys. Rev. A, vol. 136, pp. A405–A411, 1964.

[3]

I. Yarovsky and E. Evans, “Computer simulation of structure and properties of crosslinked polymers: application to epoxy resins,” Polymer (Guildf)., vol. 43, no. 3, pp. 963–969, 2002.

[5]

AC C

[4]

EP

[1]

C. Wu and W. Xu, “Atomistic molecular modelling of crosslinked epoxy resin,” Polymer (Guildf)., vol. 47, no. 16, pp. 6004–6009, Jul. 2006.

C. Wu and W. Xu, “Atomistic molecular simulations of structure and dynamics of crosslinked epoxy resin,” Polymer (Guildf)., vol. 48, no. 19, pp. 5802–5812, Sep. 2007.

[6]

J. L. Tack and D. M. Ford, “Thermodynamic and mechanical properties of epoxy resin DGEBF crosslinked with DETDA by molecular dynamics,” J. Mol. Graph. Model., vol. 26, no. 8, pp. 1269–1275, Jun. 2008.

[7]

V. Varshney, S. S. Patnaik, A. K. Roy, and B. L. Farmer, “A molecular dynamics study of epoxybased networks: Cross-linking procedure and prediction of molecular and material properties,” Macromolecules, vol. 41, no. 18, pp. 6837–6842, 2008. 17

ACCEPTED MANUSCRIPT N. Nouri and S. Ziaei-Rad, “A molecular dynamics investigation on mechanical properties of cross-linked polymer networks,” Macromolecules, vol. 44, pp. 5481–5489, 2011.

[9]

S. Yang and J. Qu, “Computing thermomechanical properties of crosslinked epoxy by molecular dynamic simulations,” Polymer (Guildf)., vol. 53, no. 21, pp. 4806–4817, Sep. 2012.

[10]

T. W. Sirk, K. S. Khare, M. Karim, J. L. Lenhart, J. W. Andzelm, G. B. McKenna, and R. Khare, “High strain rate mechanical properties of a cross-linked epoxy across the glass transition,” Polymer (Guildf)., vol. 54, no. 26, pp. 7048–7057, 2013.

[11]

M. L. Williams, R. F. Landel, and J. D. Ferry, “The Temperature Dependence of Relaxation Mechanisms in Amorphous Polymers and Other Glass-forming Liquids.,” J. Am. Chem. Soc., vol. 77, no. 14, pp. 3701–3707, 1955.

[12]

A. Soldera and N. Metatla, “Glass transition of polymers: Atomistic simulation versus experiments,” Phys. Rev. E - Stat. Nonlinear, Soft Matter Phys., vol. 74, no. 6, pp. 1–6, 2006.

[13]

A. Shokuhfar and B. Arab, “The effect of cross linking density on the mechanical properties and structure of the epoxy polymers: molecular dynamics simulation.,” J. Mol. Model., vol. 19, no. 9, pp. 3719–31, Sep. 2013.

[14]

J. D. Monk, J. B. Haskins, C. W. Bauschlicher, and J. W. Lawson, “Molecular Dynamics Simulations of Phenolic Resin: Construction of Atomistic Models,” Polymer (Guildf)., vol. 62, pp. 39–49, 2015.

[15]

N. J. Soni, P. H. Lin, and R. Khare, “Effect of cross-linker length on the thermal and volumetric properties of cross-linked epoxy networks: A molecular simulation study,” Polym. (United Kingdom), vol. 53, no. 4, pp. 1015–1019, Feb. 2012.

[16]

A. Bandyopadhyay, 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 (Guildf)., vol. 52, no. 11, pp. 2445–2452, 2011.

[17]

M. Tsige and M. J. Stevens, “Effect of crosslinker functionality on the fracture behavior of highly crosslinked polymer networks : A molecular dynamics study of epoxies,” vol. 1411, pp. 630–637, 2003.

[19]

SC

M AN U

TE D

EP

AC C

[18]

RI PT

[8]

C. Li and A. Strachan, “Molecular dynamics predictions of thermal and mechanical properties of thermoset polymer EPON862/DETDA,” Polymer (Guildf)., vol. 52, no. 13, pp. 2920–2928, 2011.

P. H. Lin and R. Khare, “Molecular simulation of cross-linked epoxy and epoxy-POSS nanocomposite,” Macromolecules, vol. 42, no. 12, pp. 4319–4327, Jun. 2009.

[20]

P. H. Lin and R. Khare, “Local Chain dynamics and dynamic heterogeneity in cross-linked epoxy in the vicinity of glass transition,” Macromolecules, vol. 43, pp. 6505–6510, 2010.

[21]

Y. Sun, Z. Zhang, K. S. Moon, and C. P. Wong, “Glass transition and relaxation behavior of epoxy nanocomposites,” J. Polym. Sci. Part B Polym. Phys., vol. 42, no. 21, pp. 3849–3858, 2004.

18

ACCEPTED MANUSCRIPT J. S. Bermejo and C. M. Ugarte, “Influence of cross-linking density on the glass transition and structure of chemically cross-linked PVA: A molecular dynamics study,” Macromol. Theory Simulations, vol. 18, no. 6, pp. 317–327, 2009.

[23]

C. K. Y. Wong, H. Fan, and M. M. F. Yuen, “Investigation of adhesion properties of Cu-EMC interface by Molecular Dynamic Simulation,” 2005.

[24]

S. Yu, S. Yang, and M. Cho, “Multi-scale modeling of cross-linked epoxy nanocomposites,” Polymer (Guildf)., vol. 50, no. 3, pp. 945–952, Jan. 2009.

[25]

M. Ionita, “Multiscale molecular modeling of SWCNTs/epoxy resin composites mechanical behaviour,” Compos. Part B Eng., vol. 43, no. 8, pp. 3491–3496, 2012.

[26]

J. Gou, B. Minaie, B. Wang, Z. Liang, and C. Zhang, “Computational and experimental study of interfacial bonding of single-walled nanotube reinforced composites,” Comput. Mater. Sci., vol. 31, no. 3–4, pp. 225–236, 2004.

[27]

J. Choi, S. Yu, S. Yang, and M. Cho, “The glass transition and thermoelastic behavior of epoxybased nanocomposites: A molecular dynamics study,” Polymer (Guildf)., vol. 52, no. 22, pp. 5197–5203, Oct. 2011.

[28]

D. Fragiadakis, P. Pissis, and L. Bokobza, “Glass transition and molecular dynamics in poly(dimethylsiloxane)/silica nanocomposites,” Polymer (Guildf)., vol. 46, no. 16, pp. 6001– 6008, 2005.

[29]

H. Sun, “Ab initio characterizations of molecular structures, conformation energies, and hydrogen-bonding properties for polyurethane hard segments,” Macromolecules, vol. 26, no. 22, pp. 5924–5936, Oct. 1993.

[30]

H. Sun, P. Ren, and J. R. Fried, “The COMPASS and validation,” vol. 8, no. l, pp. 229–246, 1998.

[31]

M. J. McQuaid, H. Sun, and D. Rigby, “Development and validation of COMPASS force field parameters for molecules with aliphatic azide chains.,” J. Comput. Chem., vol. 25, no. 1, pp. 61–71, Jan. 2004.

[32]

H. Sun, “COMPASS: An ab Initio Force-Field Optimized for Condensed-Phase Applications s Overview with Details on Alkane and Benzene Compounds,” vol. 5647, no. 98, pp. 7338– 7364, 1998.

[34]

SC

M AN U

TE D

EP

AC C

[33]

RI PT

[22]

“BIOVIA Material Studio 7.0, http://accelrys.com/products/materials-studio/.” .

A. a. Samoletov, C. P. Dettmann, and M. a J. Chaplain, “Thermostats for ‘slow’ configurational modes,” J. Stat. Phys., vol. 128, no. 6, pp. 1321–1336, 2007.

[35]

B. Leimkuhler, E. Noorizadeh, and O. Penrose, “Comparing the Efficiencies of Stochastic Isothermal Molecular Dynamics Methods,” J. Stat. Phys., vol. 143, no. 5, pp. 921–942, 2011.

[36]

H. J. C. Berendsen, J. P. M. Postma, W. F. van Gunsteren, a DiNola, and J. R. Haak, “Molecular dynamics with coupling to an external bath,” J. Chem. Phys., vol. 81, no. 8, pp. 3684–3690, 1984.

19

ACCEPTED MANUSCRIPT G. J. Martyna, D. J. Tobias, and M. L. Klein, “Constant pressure molecular dynamics algorithms,” J. Chem. Phys., vol. 101, no. 1994, p. 4177, 1994.

[38]

M. Parrinello, a. Rahman, and R. a. Parrinello M, “Crystal Structure and Pair Potentials: A Molecular_Dynamics Study,” Phys. Rev. Lett., vol. 45, no. 14, pp. 1196–1199, 1980.

[39]

T. G. Fox and P. J. Flory, “Second-order transition temperatures and related properties of polystyrene. I. Influence of molecular weight,” J. Appl. Phys., vol. 21, no. 6, pp. 581–591, 1950.

[40]

L. M. McGrath, R. S. Parnas, S. H. King, J. L. Schroeder, D. a. Fischer, and J. L. Lenhart, “Investigation of the thermal, mechanical, and fracture properties of alumina-epoxy composites,” Polymer (Guildf)., vol. 49, no. 4, pp. 999–1014, 2008.

[41]

L. Shan, K. N. E. Verghese, C. G. Robertson, and K. L. Reifsnider, “Effect of Network Structure of Epoxy DGEBA-Poly ( oxypropylene ) diamines on Tensile Behavior,” no. Epon 828, pp. 2815–2819, 1999.

[42]

B. Burton, D. Alexander, H. Klein, a Garibay-Vasquez, a Pekank, and C. Henkee, “Epoxy formulations using Jeffamine® polyetheramines,” Huntsman Corp., 2010.

[43]

A. Lee and G. B. Mckenna, “Effect of crosslink density on physical ageing of epoxy networks,” vol. 29, no. February, 1988.

[44]

D. N. Theodorou and U. W. Suter, “Atomistic modeling of mechanical properties of polymeric glasses,” Macromolecules, vol. 19, no. 1, pp. 139–154, 1986.

[45]

D. Brown and J. H. R. Clarke, “Molecular dynamics simulation of an amorphous polymer under tension. 1. Phenomenology,” Macromolecules, vol. 24, no. 8, pp. 2075–2082, 1991.

[46]

D. H. Tsai, “The virial theorem and stress calculation in molecular dynamics,” J. Chem. Phys., vol. 70, no. 3, p. 1375, 1979.

AC C

EP

TE D

M AN U

SC

RI PT

[37]

20