Molecular modeling of crosslink distribution in epoxy polymers

A small-scale MD model was established for a mixture of EPON 862 monomer (di-glycidyl ether of bisphenol-F) and crosslinking agent DETDA (diethylene toluene diamine). The molecules of EPON-862 and DETDA are shown in figure 1. A stoichiometric mixture (2 : 1 ratio) of two molecules of EPON-862 and one molecule of DETDA was modeled using the LAMMPS (large scale atomic/molecular massively parallel simulator) simulation program [17]. The initial atomic coordinates were written in a coordinate file in the native LAMMPS format and the OPLS united-atom force field [18–20] was used to define the bond, angle and dihedral parameters with Lenard-Jones non-bonded interactions. The initial 2 : 1 structure was simulated in a 10 × 10 × 10 angstrom (Å) simulation box with periodic boundary conditions. This structure was subjected to four energy minimizations (MM) and three MD simulations to drive the system to a stable equilibrium. After equilibration was achieved, the structure was replicated to form eight more structures within the simulation box such that a 16 : 8 molecular mixture of EPON 862 and DETDA monomers was established with a total of 664 united atoms. A gradual equilibration procedure was performed over a cycle of 20 MM and 10 MD simulations to equilibrate the structure and reduce the simulation box size such that a mass density of 1.2 g cm⁻³ was achieved. This procedure was previously established [16], however, it will be briefly explained here. Incremental reductions in the simulation box volume were followed by a MM simulation, a MD simulation, then another MM simulation to minimize the energy for each incremental reduction. Minimization was assumed to occur when the potential energy and pressure of the simulation boxes was stable with respect to simulated time. A total of 10 of such reductions were performed until the simulation box volume was reduced enough to yield a mass density of 1.2 g cm⁻³. Using 10 incremental volume reductions, the target mass density could be achieved without significant increases in potential energy and pressure of the molecular model. All MD simulations were conducted in the NVT (constant volume and temperature) ensemble for 100 ps at 600 K. This simulation time ensured that no residual pressures (and thus stresses) remained in the structure after each increment of simulation box size reduction. Therefore, there were no significant stresses present in the molecular structure once the density of 1.2 g cm⁻³ was achieved. The NVT ensemble made use of the Nose/Hoover thermostat for temperature control [21].

Zoom In Zoom Out Reset image size

After equilibration of the 16 : 8 molecular mixture, the structure was crosslinked based on the root mean square (rms) distance between the N atoms of DETDA and CH₂ groups of the EPON 862 molecules, similar to the approach used by Bandyopadhyay et al [16]. A crosslinked and equilibrated 16 : 8 molecular mixture is shown in figure 2 (Bonds have not been shown for better clarity of the image).

Zoom In Zoom Out Reset image size

The simulated crosslinking process is depicted in figures 3 and 4. Three assumptions were made for the crosslinking procedure:

• (1)  
  Both primary amines in DETDA were assumed to have the same reactivity.
• (2)  
  The CH₂–O and N–H bonds were broken simultaneously (figure 3).
• (3)  
  One N atom was partially activated when it had only one activated CH₂ within a defined cutoff distance.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The starting point of the crosslinking reaction is shown in figure 3 where the nucleophilic amine of the DETDA molecule reacts with the unsubstituted epoxy C, and the adjacent oxygen atom attains a negative charge when the C–O bond is broken. After forming a bond with the C atom, the N atom attains a positive charge and thus the neutrality of the EPON 862-DETDA system is maintained. Crosslinks were formed by computing all rms distances between each N atom and the CH₂ united atoms within a defined cutoff distance. The CH₂ radicals located outside the cut-off distance of a particular NH₂ group were not crosslinked to that particular group. The cutoff distance was chosen to achieve a desired level of crosslinking as described elsewhere [16]. In the following step, the H atoms were formed by breaking NH₂ bonds and were reacted with the O atoms of the broken epoxide ends. This bond formation was also performed based on the closest rms distances between the O and H atoms. The second step of the crosslinking reaction is shown in figure 4. The crosslink density of the molecular clusters is defined as the ratio of the total number of crosslinks that were formed to the 32 total crosslinks that could potentially be formed. For example, an epoxy network having 16 out of 32 crosslinks is defined as having a 50% crosslink density. A total of three molecular clusters were crosslinked to crosslink densities of 50%, 59% and 69%.

Once the molecular clusters were crosslinked and equilibrated, they were replicated to form a 3 × 3 × 3 array, such that a series of larger-scale models were established that contained 27 individual molecular cluster models. A total of three larger-scale models were initially formed, one for each of the three cluster crosslink densities (50%, 59%, and 69%). The large-scale models had 432 molecules of EPON 862 and 216 molecules of DETDA. The three large-scale systems were equilibrated in the same manner as the molecular clusters to ensure a thorough mixing of the individual molecular clusters and an overall mass density of 1.2 g cm⁻³. A total of 17 928 united atoms were used in the model, which represented a total of 25 272 real atoms. A pictorial representation of how the large-scale structure was evolved from molecular clusters is shown in figure 5.

Zoom In Zoom Out Reset image size

For each of the three large-scale molecular models, various levels of additional crosslinks were added to previously un-reacted reactive groups that spanned different molecular clusters within the same large-scale structure. Figure 6 shows the crosslinking of the large-scale molecular model after equilibration. The yellow lines in the figure are not actual crosslinked bonds but are a magnified representation in order to explain the presence of crosslinks. These additional crosslinks served to mechanically stabilize the structure and to achieve a varying degree of overall crosslink densities. These interconnecting crosslinks simulated the dispersed phase of low-crosslink density as identified in the experiments discussed above [2, 4–7]. The crosslinked clusters represented the regions of high-crosslink density and higher molecular weight. Table 1 summarizes the nine different large-scale molecular models that were established, along with their corresponding cluster crosslink densities and overall crosslink densities. The ratio of the two crosslink measures is defined as

$$\begin{equation*} R=\frac{\mbox{Cluster crosslink density}}{\mbox{Overall crosslink density}}. \end{equation*}$$

The ratio R provides a quantitative measure of the distribution of crosslinks. As shown in table 1, R varies from 0.66 to 0.99 for the large scale models. Larger values of R indicate that a higher percentage of the crosslinks are confined to the clusters, whereas small values indicate that the crosslinks are more distributed between clusters. It is also important to note that there is one 54% overall crosslinked system, two 63% overall crosslinked systems, three 70% overall crosslinked systems and three 76% overall crosslinked systems. Therefore, the predicted properties can be compared in terms of cluster crosslink density, overall crosslink density and R.

Zoom In Zoom Out Reset image size

For each of the nine cross-linked epoxy models, NPT (constant pressure and constant temperature) simulations were run for 400 ps at 29 different temperatures ranging from −73 °C (200 K) to 277 °C (550 K) at temperature intervals of 12.5 °C and at pressures of 1 atmosphere (atm). These simulations modeled the process of constant heating of the epoxy systems from cryogenic to elevated temperatures, and provided the equilibrated volumes corresponding to each temperature interval. Shown in figure 7 is the volume shrinkage of four of the large-scale models with a 50% cluster crosslink density (systems 1–4 in table 1) with respect to the volume at 277 °C. Similarly, the volume shrinkage curves of the three large-scale models with a 59% cluster crosslink density (systems 5–7 in table 1) are shown in figure 8, and the volume shrinkage curves of the two large-scale models with a 69% cluster crosslink density (systems 8 and 9 in table 1) are shown in figure 9. Data from the final 350 ps of each NPT simulation were used to establish the data points in the graphs. The data gathered during the first 50 ps of the simulations were not used because it was observed that the molecular structures were not initially equilibrated after the imposed changes in temperature, which is expected when the temperature of an initially equilibrated structure is changed. It took 50 ps for the molecular structure to reach an apparent equilibrium during the simulations. This time interval was determined by observing the oscillation of the total pressure. The oscillation arrived at a constant level by 50 ps. Figures 7–9 also show the scatter present in the final 350 ps of each simulation in the form of error bars. It is observed that the scatter is relatively small at low temperatures and large at high temperatures, which indicates that the crosslinked epoxy polymer molecules undergo larger vibrations in the rubbery phase to cause the changes commonly associated with the T_g.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The data in figures 7–9 show that for decreasing simulation temperatures, the volume shrinkage increases, as expected, due to thermal contraction. The data also clearly show that for increasing levels of overall crosslink density the volume shrinkage decreases. This trend holds for all three cluster crosslink densities, and is similar to that reported by Bandyopadhyay et al [16].

The epoxy systems shown in figure 7 differ from each other in one aspect; the number of crosslinks in between the crosslinked clusters. As these crosslinks connecting the clusters are considered as the dispersed phase of low crosslink density, figure 7 shows that the clusters become more rigidly bonded to each other as the dispersed phase becomes more crosslinked. Thus the entire system becomes more resistant to volumetric deformations for a given temperature change. These results are expected, and are presented here as a validation of the present modeling technique.

Linear regression lines were fitted on the first 10 data points (−73 to 39.5 °C) and the last 10 data points (164.5 to 277 °C) of the volume shrinkage curves shown in figures 7, 8 and 9 for determining the coefficients of volumetric thermal expansion (CVTE) below and above T_g, respectively. The T_g of similar EPON 862-DETDA systems have been reported elsewhere [12, 14, 16, 22]. The remaining nine data points in the middle portion of the volume shrinkage curves were assumed to be the glass transition range and were not used in calculating CVTE. The CVTE is defined as

$$\begin{equation} \beta =\frac{1}{V_0 }\left( {\frac{\partial V}{\partial T}} \right)_P, \end{equation} \tag{ 1 }$$

where V₀ is the initial volume of the simulation box at 277 °C, and the subscript P implies a constant-pressure process. The coefficient of linear thermal expansion (CLTE) is defined as

$$\begin{equation} \alpha =\frac{\beta }{3}=\frac{1}{L_0 }\left( {\frac{\partial L}{\partial T}} \right)_P , \end{equation} \tag{ 2 }$$

where L₀ is the initial length of each side of the cubic simulation box at 277 °C. The CLTE values obtained for the nine different crosslinked systems are shown in figure 10. The uncertainty in the predicted CLTE values is shown in figure 10 based on the uncertainty in the data shown in figures 7, 8, and 9 and standard rules of propagation of errors [23].

Zoom In Zoom Out Reset image size

The predictions shown in figure 10 are consistent with those described in the literature [12, 16, 24], that is, the CLTE values decrease with increasing levels of overall crosslink density, both above and below T_g. This is likely because of increasing numbers of covalent crosslinks in the polymer network serve to inhibit volume expansion in response to temperature increases. The larger values of CLTE above T_g are due to the increase in molecular mobility of the polymer network in the rubbery state. In the glassy state, the network has a higher density and thus reduced molecular mobility decreases the material's volumetric response to temperature changes.

The CLTE values of the systems with overall crosslink densities of 63%, 70% and 76% are shown in tables 2, 3, and 4, respectively, with the corresponding values in uncertainty, calculated as described above. The arrangement of the data in this manner allows for the direct comparison between systems with the same overall crosslink density and different crosslink distributions. The comparison of systems with the same overall crosslink density but differently crosslinked clusters provides insight into the effect of crosslink distribution. Each of the three tables shows variation in values of CLTE above and below T_g. The scatter shown in figures 7–9 had a small effect on the CLTE values as seen in tables 2–4.

The CLTE data are plotted against R in figure 11. Above T_g, there is a considerable amount of scatter in the data, and there is no clear influence of R on the CLTE. Below T_g, the scatter in the data is reduced and it appears that there may be a slight increase in the CLTE for increases in R. This would indicate that as more crosslinks are confined to the individual clusters, the overall molecular mobility of the system increases by a small amount, as there is less resistance to volumetric changes with respect to temperature changes.

Zoom In Zoom Out Reset image size

Because the data in figure 11 indicate there is a slight influence of crosslink distribution on the overall CLTE (no influence above T_g), it follows that the scatter in the data shown in figure 10 is not highly influenced by the crosslink distribution, especially for the above-T_g data. Therefore, it is likely that the scatter in both sets of data in figure 10 is mostly caused from a statistical variation in properties of the discrete systems. Each system has a unique molecular structure and a finite number of atoms, and thus it is expected that a finite amount of scatter will exist in predicted properties. This has been discussed in detail elsewhere for thermoplastic systems [25, 26]. It is expected that as the number of atoms increases, the scatter will decrease.

The elastic constants of the nine crosslinked systems were predicted with MD simulations. Various magnitudes of uniaxial tensile strain and volumetric strain were applied to the MD models to determine the elastic properties. All simulations were performed in the NVT ensemble under periodic boundary conditions. For all simulations, the applied strains were ramped up to the final strain level over the duration of the simulated time. The virial stress components were calculated for each timestep over the entire loading range. It was assumed that the bulk material mechanical response was linear-elastic and isotropic. It is important to note that although the elastic response of polymers is generally dependent on applied strain rate, and the simulated deformations are at relatively high strain rates, it is expected that the predicted elastic properties are close to those observed experimentally [15, 16]. This has been demonstrated for a range of crosslink densities [16]. Therefore, it is expected that the predicted mechanical response from the high- and low-crosslink density regions are not dependent on the simulated strain rate.

For the volumetric deformations, equal axial strains were applied simultaneously over all three principal directions of the simulation boxes using the 'fix deform' and 'fix nvt/sllod' commands provided by the LAMMPS software. The simulations were run at 300 K (room temperature) for 200 ps with timesteps of 0.2 fs. The simulated deformations were applied in both tension and compression according to

$$\begin{equation} \varepsilon _{xx} =\varepsilon _{yy} =\varepsilon _{zz} =\pm 0.003,\pm 0.005,\pm 0.007 \end{equation} \tag{ 3 }$$

where ε_xx, ε_yy and ε_zz are the axial components of the infinitesimal strain tensor component with respect to the basis vectors defined in figure 5. The overall dilatation of the molecular models was

$$\begin{equation} \Delta =\varepsilon _{xx} +\varepsilon _{yy} +\varepsilon _{zz} . \end{equation} \tag{ 4 }$$

For each timestep in these simulations, the overall hydrostatic stress σ_h of the model was calculated by

$$\begin{equation} \sigma _{\rm h} =\textstyle{1 \over 3}(\sigma _{xx} +\sigma _{yy} +\sigma _{zz} ), \end{equation} \tag{ 5 }$$

where σ_xx, σ_yy, and σ_zz are the volume-averaged virial stress tensor components. For each timestep of the final 180 ps of the simulations, the hydrostatic stress and dilatation were used to perform a linear regression analysis to determine the bulk modulus

$$\begin{equation} K=\frac{\sigma _{\rm h} }{\Delta }. \end{equation} \tag{ 6 }$$

The bulk moduli calculated for positive and negative dilatations over each of the three strain levels were averaged. The data gathered during the first 20 ps of the simulations were not used in the regression analysis because it was observed that the molecular structures were not initially equilibrated after the imposed deformation, which is expected when an initially equilibrated structure is statically perturbed. It took 20 ps for the molecular structure to reach an apparent equilibrium during the simulated deformation. This time interval was determined by observing the oscillation of the total pressure. The oscillation arrived at a constant level by 20 ps. A total of nearly one million data points were used to calculate the bulk modulus for each epoxy system.

For the uniaxial tensile deformations, strains were individually applied along each of the principal axes of the simulation box using the 'fix deform' and 'fix nvt/sllod' commands in the LAMMPS software package. Axial strains were applied such that Poisson contractions were allowed in the transverse directions. For each model, the tensile strain was divided equally over 200 ps with timesteps of 0.4 fs. Tensile strains of magnitudes ±0.003, ±0.005 and ±0.007 were applied in each direction.

A linear regression analysis of stress versus strain was performed over the final 180 ps of each simulation to determine the corresponding Young's modulus (E). Young's moduli were averaged over all three axes and all three strain levels in both tension and compression. Using the values obtained for Bulk moduli and Young's moduli, Shear moduli (G) were calculated for all nine molecular systems as

$$\begin{equation} G=\frac{3KE}{9K-E}. \end{equation} \tag{ 7 }$$

The predicted Young's and shear moduli are shown in figure 12 for the entire range of overall crosslink densities. From the figure it is clear that there is an increasing trend in both moduli as the crosslink density increases. The magnitude of the predicted elastic properties is similar to those reported in the literature for the same epoxy system [15, 16, 27–29]. This trend with respect to overall crosslink density is similar to that observed by Bandyopadhyay et al [16] and Li and Strachan [15] for EPON 862/DETDA. Similar increasing trends have been reported in other epoxy systems by Gupta et al [27] and in epoxy–nanotube composites by Miyagawa et al [29]. Lees and Davidson [28] discuss data that show increasing distances between reactive sites in epoxies result in decreasing crosslink densities and decreasing elastic properties. Figure 12 also indicates that there is a substantial amount of scatter in the predicted elastic properties over the range of overall crosslink densities, particularly at overall crosslink densities of 63% and 70%.

Zoom In Zoom Out Reset image size

Tables 5, 6 and 7 show Young's modulus and shear modulus of the systems with overall crosslink densities of 63%, 70% and 76%, respectively. Similar to the CLTE values presented in tables 2, 3 and 4, there is no clear trend in the elastic properties with changes in the crossink distribution. The lack of trend may be due to the significant amount of scatter in the data, which is substantially reduced for the system with the overall crosslink density of 76%.

Figure 13 shows Young's and shear moduli plotted as a function of R. The regression lines show R has no influence on moduli values. Similar to the interpretation of figure 11, this indicates that the crosslink distribution has no detectable influence on the overall elastic response. The scatter in the data shown in figures 12 and 13 is therefore most likely due to statistical variation in properties of the discrete molecular systems. This scatter would likely disappear if it were possible to simulate molecular systems with many more orders of magnitude of atoms.

Zoom In Zoom Out Reset image size