# Molecular Simulations of Thermo-Mechanical Properties of Epoxy-Amine Resins

Mathilde Orselly, Julien Devemy, Agathe Bouvet-Marchand, Alain Dequidt, Cédric Loubat, Patrice Malfreyt

► To cite this version:

Mathilde Orselly, Julien Devemy, Agathe Bouvet-Marchand, Alain Dequidt, Cédric Loubat, et al. Molecular Simulations of Thermo-Mechanical Properties of Epoxy-Amine Resins. ACS Omega, 2022, 10.1021/acsomega.2c03071. [PDF] [ DOI ]

HAL Id: hal-03759439

https://uca.hal.science/hal-03759439v1

Submitted on 24 Aug 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
<figure><img src="image_2.png" /></figure>

# Molecular Simulations of Thermomechanical Properties of Epoxy-Amine Resins

Mathilde Orselly,* Julien Devemy, Agathe Bouvet-Marchand, Alain Dequidt, Cédric Loubat, and Patrice Malfreyt*

<figure><img src="image_3.png" /></figure>

Cite This: ACS Omega 2022, 7, 30040–30050

<figure><img src="image_4.png" /></figure>

Read Online

## ACCESS

<figure><img src="image_5.png" /></figure>

Metrics & More

<figure><img src="image_6.png" /></figure>

Article Recommendations

<figure><img src="image_7.png" /></figure>

Supporting Information

**ABSTRACT:** All-atom molecular dynamics (MD) simulations were performed with the CHARMM force field to characterize various epoxy resins, such as aliphatic and bisphenol-based resins. A multistep cross-linking algorithm was established, and key properties such as density, glass temperature, and elastic modulus were calculated. A quantitative comparison was made and was proven to be in good agreement with experimental data, with average absolute deviations between experiments and molecular simulation comprised between 2% and 12%. Additional findings on structure–property relationships were highlighted such as the effect of the cross-linking rate and oligomerization of the resin.

<figure><img src="image_8.png" /></figure>

## 1. INTRODUCTION

Epoxies are some of the most prominent thermosetting polymers valued nowadays due to their wide range of applications. This class of polymers is of interest for instance for structural components in aerospace applications, as well as electronics packaging or various types of coatings. This is made possible by their interesting versatility. The final material can be easily modulated by the chemical structure of the epoxy or the nature of the hardener, responsible for the striking thermal and mechanical properties of the resulting material.¹

Therefore, it is of great importance to be able to characterize them and comprehend the influence of the structural network on their properties. Aside from the usual experimental testing, computational chemistry has proven to be a great asset in this task. Indeed, predictive simulations are playing an increasing role in understanding the atomistic origin of the properties of polymers as well as in speeding up the process of selecting materials, thus complementing experiments in the search for tailored materials. However, one of the greatest challenges is to be able to reproduce reliably the physical, mechanical, and thermodynamic properties of such materials from atomistic models. Molecular dynamics (MD) has proven to be a powerful method to generate equilibrated atomistic configurations on which macroscopic properties are averaged.²

One of the first steps to reproduce thermosetting polymers in molecular simulations is to be able to cross-link the monomers, making it possible to characterize the polymeric material. Two approaches are typically used for the cross-linking of the monomers: a “back-mapping” procedure$^{3–6}$ and a “reactive” one.$^{7–15}$ The former is based on the cross-linking of coarse-grained monomers and then the all-atom model is

later back-mapped in lieu of the coarse grains, while the latter already considers an all-atom model and the cross-linking is done between the experimental reactive pairs. These cross-linking methods can also be single or multiple steps but exhibit the same topology overall.¹⁶

Regarding the epoxy systems investigated by molecular simulations, the most extensive research effort was found in resins resulting from the cross-linking between the diglycidyl ether of bisphenol F (DGEBF) and diethyl toluenediamine (DETDA),$^{9,11,12,17,18}$ triethylenetetramine (TETA)$^{10,19–21}$ or diethylenetriamine (DETA),$^{22}$ especially in the fields of modern aeronautics. The cross-linking of the diglycidyl ether of bisphenol A (DGEBA) resin was also widely studied with different hardeners, such as isophorone diamine (IPDA),$^{23–25}$ trimethylene glycol di-*p*-aminobenzoate (TMAB),$^{26}$ diethyltoluenediamine (DETDA),$^{4,27–29}$ triethylenetetramine (TETA),$^{27}$ ethylenediamine (EDA),$^{7}$ diaminodiphenyl sulfone (DDS),$^{30}$ methylenedianiline (MDA),$^{31}$ poly(oxopropylene) diamines (POP),$^{16,32,33}$ 4,4'-methylenebis(cyclohexylamine) (MCA),$^{33}$ diethylenetriamine (DETA),$^{13,34–36}$ and polyether-amine JEFFAMINE D-230.$^{37–40}$

To the best of our knowledge, novel epoxy resins based on the diglycidyl ether of vanillyl alcohol (DGEVA), diglycidyl ether of 1,4-butanediol (DGEBU), triglycidyl ether of

Received: May 17, 2022

Accepted: August 1, 2022

Published: August 22, 2022

<figure><img src="image_9.png" /></figure>
<figure><img src="image_10.png" /><figcaption>**Figure 1.** Molecular representation of the epoxy resins (a–e) and of the hardeners (f–m). Carbon, oxygen, nitrogen, and hydrogen atoms are shown in gray, red, blue, and white, respectively. (n, o) Snapshots of cross-linked systems, where orange and gray phases represent respectively the epoxy resins and hardeners.</figcaption></figure>

phloroglucinol (PHTE), 4,4'-methylenebis(2,6-diethylaniline) (MDEA), *m*-xylylenediamine (MXDA), 4,4'-methylenebis(2-isopropyl-6-methylaniline) (MMIPA), tetraethylenepentamine (TEPA), and 4,7,10-Trioxa-1,13-tridecane diamine (TTDA) have not yet been studied by molecular simulations. Due to the broad variety of chemical structures, a comparison between experimental and simulation results on several key properties such as density, glass transition temperature, and elastic moduli is made possible.

Through this, an evaluation of the performance and the reliability of general atomistic force fields and the capacity of simulation methods to reproduce the thermomechanical properties of these resins is carried out. Indeed, predicting the physical, mechanical, and thermodynamic properties of resins by using full atomistic models remains challenging because of the complexity of the interactions involved and the algorithm needed to post-treat the data. Such molecular simulations can open the way to interpreting the thermomechanical properties of resins on the basis of their structure and the molecular interactions involved.

In order to establish this comparison, a cross-linking algorithm was first made based on previous studies using the reactive method.⁷⁻¹⁵ This was followed by several procedures to ensure the reliability of the properties obtained, thus making it possible to enlighten key properties of materials and their relationship to their structure. To some extent, this study also aims at evaluating the transferability of the molecular models on resins not yet investigated by molecular simulations.

The paper is organized as follows. **Section 2** describes the experimental and computational procedures. **Section 3** presents the main results of this work in terms of density, glass transition temperature, and Young's modulus of various epoxy-amine resins. We conclude in **Section 4**.

## 2. EXPERIMENTAL AND SIMULATION METHODS

### 2.1. Experimental Details. 2.1.1. Hardeners and Resins Studied.
Various epoxy resins are studied, exhibiting a large panel of thermomechanical properties. These resins are widely used in industry. The formulations presented in this study include the compounds DGEBA, DGEBF, DGEBU, PHTE, and DGEVA as well as the hardeners IPDA, DETA, EDA, TEPA, TTDA, MDEA, MMIPA, and MXDA. DGEVA and PHTE can be derived from bio-based raw materials and were synthesized by Specific Polymers. These components are presented in Figure 1. The average number of repeating units (n) of DGEBA (Sika, Biresin CR170) was determined by ¹H NMR (Figure 2).

<figure><img src="image_11.png" /><figcaption>**Figure 2.** Depiction of the meaning of *n* in terms of the oligomerization of DGEBA.</figcaption></figure>

**Table 1** gives for the epoxies studied the degree (n) of oligomerization of the resin, the resin:hardener ratio, the molecular weight of the constitutive repeat unit, and the empirical density (detailed in the following subsection) as well as the number of atoms present and the size of the box for the simulations.

### 2.1.2. Determination of the Experimental Density and Glass Transition Temperature.
The measurement of density has been performed on materials fully cured using a Ohaus density determination kit for a scale. Density determination is performed by Archimedes' principle. The density of the material is determined with the aid of a liquid whose density $ρ_0$
<table><caption>Table 1. Description of the Characteristics of Resins and Hardeners Studied, Breaking down the Chemical Structure of the Repeat Unit: Composition, Mass, and Expected Density<sup>a</sup></caption><thead><tr><th>resin</th><th>hardener</th><th>n</th><th>R:H ratio</th><th>mass of CRU (g/mol)</th><th>density (g/cm<sup>3</sup>) (eq2)</th><th>no. of atoms</th><th>side of the box (Å)</th></tr></thead><tbody><tr><td>DGEBU</td><td>IPDA</td><td>0.32</td><td>2:1</td><td>574</td><td>1.040</td><td>7938</td><td>41.9</td></tr><tr><td>DGEBA</td><td>DETA</td><td>0.32</td><td>5:2</td><td>2368</td><td>1.141</td><td>11267</td><td>48.7</td></tr><tr><td>DGEBA</td><td>EDA</td><td>0.32</td><td>2:1</td><td>924</td><td>1.145</td><td>10919</td><td>47.8</td></tr><tr><td>DGEBA</td><td>IPDA</td><td>0.32</td><td>2:1</td><td>1034</td><td>1.117</td><td>12701</td><td>50.5</td></tr><tr><td>DGEBA</td><td>MMIPA</td><td>0.32</td><td>2:1</td><td>1174</td><td>1.128</td><td>14240</td><td>52.6</td></tr><tr><td>DGEBA</td><td>MXDA</td><td>0.32</td><td>2:1</td><td>1000</td><td>1.152</td><td>11729</td><td>49.1</td></tr><tr><td>DGEBA</td><td>TEPA</td><td>0.32</td><td>7:2</td><td>3404</td><td>1.137</td><td>11639</td><td>49.2</td></tr><tr><td>DGEBA</td><td>TTDA</td><td>0.32</td><td>2:1</td><td>1084</td><td>1.129</td><td>13106</td><td>50.8</td></tr><tr><td>DGEBF</td><td>IPDA</td><td>0</td><td>2:1</td><td>1008</td><td>1.140</td><td>12150</td><td>49.8</td></tr><tr><td>DGEVA</td><td>DETA</td><td>0</td><td>5:2</td><td>1616</td><td>1.174</td><td>11500</td><td>48.5</td></tr><tr><td>DGEVA</td><td>IPDA</td><td>0</td><td>2:1</td><td>734</td><td>1.134</td><td>8910</td><td>44.1</td></tr><tr><td>DGEVA</td><td>MDEA</td><td>0</td><td>2:1</td><td>874</td><td>1.146</td><td>12900</td><td>50.3</td></tr><tr><td>DGEVA</td><td>MXDA</td><td>0</td><td>2:1</td><td>700</td><td>1.187</td><td>9800</td><td>45.7</td></tr><tr><td>PHTE</td><td>IPDA</td><td>0</td><td>3:4</td><td>1686</td><td>1.119</td><td>12900</td><td>50.3</td></tr><tr><td>PHTE</td><td>MDEA</td><td>0</td><td>3:4</td><td>2106</td><td>1.136</td><td>15750</td><td>54.1</td></tr><tr><td>PHTE</td><td>MXDA</td><td>0</td><td>3:4</td><td>1584</td><td>1.186</td><td>11100</td><td>48.0</td></tr></tbody></table>

is known (water was used as an auxiliary liquid here). The solid is weighed in air (A) and then in the auxiliary liquid (B). The density $\rho$ can be calculated from the two weights as follows $\rho = [A/(A - B)](\rho_0 - \rho_L) + \rho_L$ with $\rho_0$ being the density of the auxiliary liquid and $\rho_L$ the air density (0.0012 g/cm$^3$). The temperature of the auxiliary liquid was measured and taken into account in all density determinations.

DSC analyses were performed by using a Mettler Toledo DSC instrument. Samples having masses of approximately 10 mg were placed in 100 µL aluminum pans with pierced lids under a nitrogen atmosphere. The applied heating rate was 10 °C min$^{-1}$ under a nitrogen atmosphere. The thermal behavior of the samples was investigated using two repeated heating-cooling cycles. The following thermal procedure was used: first ramp from -60 to 250 °C and then from 250 to 25 °C. The glass transition temperature ($T_g$) was determined from the second heating curve.

## 2.2. Computational Details.
The CHARMM general force field (CgenFF)$^{41-43}$ is used to describe the intra- and intermolecular interactions of the thermosetting polymers. This class I force field has the following functional form:

$$ U_{\text{intra}} = \sum_{\text{bonds}} k_r(r_{ij} - r_0)^2 + \sum_{\text{angles}} k_{\theta}(\theta_{ijk} - \theta_0)^2 \\ + \sum_{\text{dihedrals}} \sum_{n} k_{\phi,n}[\cos(n\phi_{ijkl} + \delta_n) + 1] \\ + \sum_{\text{Impropers}} k_{\phi}(\phi_0 - \phi_0)^2 + \sum_{\text{Urey-Bradley}} k_{\text{ub}}(r_{1,3} - r_{1,3_0})^2 + \sum_{\text{Coulomb}} \frac{q_i q_j}{4\pi\epsilon_0 r_{ij}} + \sum_{\text{VdW}} \epsilon_{ij} \left[ \left( \frac{R_{\text{min},ij}}{r_{ij}} \right)^{12} - 2 \left( \frac{R_{\text{min},ij}}{r_{ij}} \right)^6 \right] \quad (1) $$

For completeness, it is worth noting that force field parameters are usually denoted using $\sigma$ rather than $R_{min}$, the relationship between the two being $\sigma = \frac{R_{min}}{2^{1/6}}$. Additionally, the CHARMM force field presents an additional component in the form of the

Urey–Bradley term, which consists of a harmonic potential as a function of the distance between the nonbonded atoms 1 and 3 of a 1–2–3 angle.

The Lennard–Jones parameters between pairs of different atoms are obtained from the Lorentz–Berthelot combination rules, in which $\epsilon_{ij}$ values are calculated by using the geometric mean and $R_{min}$ values by using the arithmetic mean. The cutoff distance for the van der Waals interactions is set to 8.0 Å, while the particle–particle particle mesh (PPPM) method with an accuracy of $1.0 \times 10^{-4}$ and a cutoff of 8 Å is considered for the long-range electrostatic interactions. For both the Lennard–Jones and Coulombic interactions, a force-switching function was applied between 8 and 10 Å.

### 2.2.1. Establishing the Cross-Linking Algorithm.
In the scope of this work, a multistep reactive cross-linking algorithm was used and is described in Figure 3, where the statistical NVT and NPT ensembles refer to the canonical and isothermal–isobaric ensembles, respectively. It includes the following steps, where the sequence of 3 + 4 + 5 steps was repeated until a predefined conversion rate was acquired (here 95% in accordance with the experimental conversion obtained after curing).

1. Stoichiometric packing of the unreacted monomers into a box of predefined density using Packmol.$^{44}$ The parameters (type, bond, angle, dihedral, charges) are those of the reacted mixture.

2. Relaxation of the initial mixture with an NPT ensemble during 500 ps at 300 K.

3. NPT simulation at 300 K to relax the system.

4. NVT simulation at 900 K to allow the quick movement of the atoms while searching for unreacted pairs within a defined cutoff distance. At each time step, the interatomic distances between unreacted pairs are calculated. The simulation was stopped once a pair was found within 2 Å.

5. Creation of the bonds for two pairs that are in agreement with the distance conditions. Additionally, a second cutoff of 3 Å was applied to allow the formation of multiple bonds per loop and to quicken the process.

6. Updating of the parameters for the pairs that did not react.

<footnote><sup>a</sup>Calculated from eq 2 as the number of atoms present and the size of the box for the simulations.</footnote>
<figure><img src="image_12.png" /><figcaption>**Figure 3.** Multistep reactive cross-linking algorithm.</figcaption></figure>

7. Quick NPT relaxation at 300 K on the fully cross-linked system during 1 ns.

8. NPT simulation during 2 ns at 600 K to further relax the system, following a ramp from 300 to 600 K for 1 ns.

9. Last NPT simulation during 2 ns at 300 K to obtain a fully relaxed system, following a ramp from 600 to 300 K for 1 ns.

This cross-linking algorithm was designed without regard to the actual polymerization reaction kinetics. Experimentally, the curing temperature of such systems is below 200 °C. However, here the cross-linking was done partially at 627 °C (above the thermal decomposition temperature). The need for a higher temperature was driven by the appeal of a quick rearrangement of the atoms within the span of the MD simulations. Additionally, a buck pair style was used for the unreacted pairs. Only the attractive term of the Lennard-Jones potential ($1/r^6$) was used. Reactive pairs are thus attracted to one another in order to facilitate bond creations. Without this potential between reactive pairs, the acquisition of a highly cross-linked resin as depicted in this work would be very time-consuming. Additionally, the cutoff range was chosen in agreement with Figure 4. Radial distribution functions, shown in Figure 4, were calculated on an uncured mixture of DGEBA +IPDA and DGEVA+MXDA between the reactive carbon of the epoxy resin and the reactive nitrogen of the amine hardener. The first maxima are close to 3.7 and 4.1 Å, respectively. This is why a cutoff of 2–3 Å is deemed relevant when it comes to the compromise between computational efficiency and the resulting cross-linked resins. Too high of a cutoff guarantees faster cross-linking but can also lead to the

<figure><img src="image_13.png" /><figcaption>**Figure 4.** Radial distribution functions of uncured resin DGEBA +IPDA (red) and DGEBA+MXDA (green) between the reactive carbon of the epoxy resin and the reactive nitrogen of the amine hardener.</figcaption></figure>

formation of metastable resins where bonds cannot come back to their equilibrium value (~1.45 Å).

Experimentally, primary amines go through an addition reaction with the epoxide groups to form by ring-opening a hydroxyl group and a secondary amine. The secondary amine can further react with an epoxide to form a tertiary amine and an additional hydroxyl group. A hydroxyl-epoxide addition can also occur but less frequently. This last step was ignored during the simulated cross-linking, as well as the difference in reactivity of the primary and secondary amines.

All potential reactive sites were activated equally by removing the hydrogen atoms from the amine of the curing agent and hydrating the oxygen atoms from the epoxide and thus opening up the oxirane. The force field parameters (bonds, angles, dihedrals, impropers, partial charges) were thus adjusted to take into account this reacted structure. This change was revoked for the molecules that did not react at the end of the cross-linking procedure (see step 6 in the list).

The sequence 7 + 8 + 9 in the list was the postcure annealing process performed to fully relax the system. Using a high-temperature annealing period after cross-linking increases the density due to the fact that this addition of thermal energy can disrupt metastable configurations and allow finding lower energy configurations closer to global minima. It was successfully demonstrated that the thermal agitation produced higher density configurations, even though the original densities had already plateaued with the first NPT equilibration.²²

Eight simulation boxes were hence created for each system containing the reacted form of the liquid mixture. They differ in their structure and the initial positioning of the molecules as to produce a statistical sample. The multistep cross-linking procedure depicted in Figure 3 was used to polymerize the resins at a conversion rate of 95% of the reactive end groups. The use of a high temperature of 900 K and a cutoff value of 2–3 Å grants a high rate of conversion while staying true to the interatomic distances. The full procedure can last between 6 and 20 ns.

To conclude on this methodological aspect, this method was employed to achieve the same degree of conversion as found experimentally within the time-limit constraints of the MD simulations. The accuracy of the structural, volumetric, and
thermal properties of the resulting cross-linked polymer will be evaluated in the course of this paper.

### 2.2.2. Description of the MD Simulations.
The velocity-Verlet integrator was used to integrate the equations of motion using a time step of 1 fs. The SHAKE algorithm was applied to constrain the length of covalent bonds to hydrogen atoms to their equilibrium values. In order to control the temperature and pressure for NVT and NPT simulations, a Nosé–Hoover thermostat and barostat were applied.⁴⁵ Unless stated otherwise, the MD simulations were performed at 300 K and 1 atm. The system was equilibrated during the first nanosecond, and data were acquired during the following nanosecond. Eight simulation boxes were used for each system. The molecular dynamics calculations were performed using the LAMMPS package.⁴⁶

## 3. RESULTS AND DISCUSSION

### 3.1. Measuring the Density.
Before examining the density of cured resins, the uncured resins and hardeners were examined to appreciate the relevance of the CHARMM force field. Figure 5 shows that the majority of the monomers are included in the 5% domain, with a warning concerning the hardener EDA.

<figure><img src="image_14.png" /><figcaption>Figure 5. Comparison of experimental and simulated densities of the uncured resins and hardeners. The dotted red lines correspond to a deviation of 5% from the experimental value.</figcaption></figure>

Figure 6 reports the values of density of different epoxy-amine resins obtained from experimental measurements, empirical calculations, and atomistic simulations. Results taken from other works are given in Table 2. An empirical power law was taken into account to provide a qualitative comparison for resins where the experimental density was not present at room temperature. This relation can be used to describe the density of linear and network polymers according to their average atomic mass ($M_a$) and constitutive repeat unit (CRU) as stated in eq 2$^{53}$

$$ \rho = k \times M_a^{2/3} \quad (\text{kg/m}^3) \qquad (2) $$

where $k \approx 320$, $\rho$ is in kg/m$^3$, and $M_a$ is in g/mol. $M_a$ is the average atomic mass of the polymer. This empirical relationship remains valid for most organic polymers ($1.1 < \rho < 1.4$ g cm$^{-3}$). Unfortunately, the maximum error of up to 10% endangers a quantitative comparison, thus requiring the

<figure><img src="image_15.png" /><figcaption>Figure 6. Comparison of the simulated, empirical, and experimental densities. Results taken from other works are given in Table 2.</figcaption></figure>

knowledge of the experimental density for a more refined comparison.

First, we observe that the agreement between empirical calculations (eq 2) and molecular simulation is good, with an average absolute deviation of 2.2%. The predicted values by molecular simulations also reproduce the variations in densities among the various resins. Other molecular models have also been shown to reproduce accurately the density for DGEBA/IPDA$^{24,25}$ and DGEBA/DETA.$^{34}$ We can conclude that the CHARMM force field is able to yield simulated densities of epoxy-amine resins that match reasonably with experiments, even if this model is not exclusively designed for polymer and additives, unlike the COMPASS model. The good agreement between the different density values obtained by simulation and empirical calculations enables a validation of the choice of the force field and the methodology used to polymerize the cross-linked epoxy resins. The CHARMM force field can contribute for epoxy resins to the determination of a qualitative study but offers an error rate of up to 5.7% with the few experimental points. The simulated values reproduce well the variations between resins but give an underestimation of the density. On the other hand, the experimental density calculated with the Archimedes method can also present an underestimation of the real density of the polymer. Taking into account this bias, this first comparison of the density shows that the model and the methodology developed for the cross-linked epoxy resins perform well for this thermodynamic property. To go further in the validation of the methodology, we propose to extend the comparison to the glass transition temperature.

### 3.2. Assessing the Glass Transition Temperature.

Modeling of glass transition temperature is based on the dilatometry principle; the density decreases with an increase in the temperature. The glass transition temperature is then characterized by a break in this slope, showing a change in the relaxation speed of the material. It separates the vitreous from the rubbery states. The modeling procedure consists of heating the polymer at $T_g+200$ K for 2 ns and cooling it at a rate of 20 K ns$^{-1}$ down to $T_g-200$ K. This cooling is done stepwise to retrieve the density of the equilibrated system at each plateau of temperature. The value of the ramp is a good compromise
<table><caption>Table 2. Compiled Simulated and Experimental Data Sets for Cross-Linked Resins at 95% Confidence<sup>a</sup></caption><thead><tr><th rowspan="2">system</th><th colspan="2">density (g/cm<sup>3</sup>)</th><th colspan="2">T<sub>g</sub> (°C)</th><th colspan="2">E (GPa)</th></tr><tr><th>exptl</th><th>sim</th><th>exptl</th><th>sim</th><th>exptl</th><th>sim</th></tr></thead><tbody><tr><td>DGEBU/IPDA</td><td></td><td>1.050<sub>3</sub></td><td>60<sup>47</sup></td><td>55<sub>8</sub></td><td></td><td>2.23<sub>11</sub></td></tr><tr><td>DGEBA/DETA</td><td>1.178</td><td>1.111<sub>3</sub></td><td>119,<sup>48</sup> 132,<sup>49</sup> 142*</td><td>123,<sup>34</sup> 133<sub>10</sub></td><td>2.55<sub>11</sub><sup>48</sup></td><td>2.57<sub>14</sub></td></tr><tr><td>DGEBA/EDA</td><td></td><td>1.125<sub>3</sub></td><td>70,<sup>50</sup> 118,<sup>7</sup> 115<sup>51</sup></td><td>62<sub>3</sub>, 117<sup>7</sup></td><td></td><td>2.79<sub>15</sub></td></tr><tr><td>DGEBA/IPDA</td><td>1.135</td><td>1.084<sub>5</sub></td><td>150*, 155,<sup>48</sup> 159<sup>49</sup></td><td>149<sub>21</sub>, 158<sup>25</sup></td><td>2.65<sub>5</sub><sup>48</sup></td><td>2.56<sub>13</sub></td></tr><tr><td>DGEBA/MMIPA</td><td>1.118</td><td>1.075<sub>3</sub></td><td></td><td>136<sub>7</sub></td><td></td><td>2.19<sub>20</sub></td></tr><tr><td>DGEBA/MXDA</td><td>1.185</td><td>1.124<sub>3</sub></td><td></td><td>101<sub>10</sub></td><td></td><td>3.41<sub>24</sub></td></tr><tr><td>DGEBA/TEPA</td><td></td><td>1.102<sub>2</sub></td><td>126<sup>48</sup></td><td>125<sub>24</sub></td><td>2.34<sub>13</sub><sup>48</sup></td><td>2.43<sub>23</sub></td></tr><tr><td>DGEBA/TTDA</td><td></td><td>1.116<sub>2</sub></td><td>69<sup>47</sup></td><td>84<sub>8</sub></td><td></td><td>2.14<sub>6</sub></td></tr><tr><td>DGEBF/IPDA</td><td></td><td>1.102<sub>4</sub></td><td></td><td>129<sub>10</sub></td><td></td><td>1.92<sub>16</sub></td></tr><tr><td>DGEVA/DETA</td><td></td><td>1.179<sub>2</sub></td><td>79*</td><td>32<sub>10</sub></td><td></td><td>2.57<sub>13</sub></td></tr><tr><td>DGEVA/IPDA</td><td></td><td>1.154<sub>6</sub></td><td>97,<sup>52</sup> 107*</td><td>89<sub>8</sub></td><td></td><td>3.99<sub>26</sub></td></tr><tr><td>DGEVA/MDEA</td><td></td><td>1.144<sub>4</sub></td><td>104*</td><td>98<sub>13</sub></td><td></td><td>3.23<sub>18</sub></td></tr><tr><td>DGEVA/MXDA</td><td></td><td>1.219<sub>3</sub></td><td>74*</td><td>75<sub>9</sub></td><td></td><td>3.41<sub>24</sub></td></tr><tr><td>PHTE/IPDA</td><td></td><td>1.108<sub>5</sub></td><td>225*, #</td><td>226<sub>15</sub></td><td></td><td>3.99<sub>33</sub></td></tr><tr><td>PHTE/MDEA</td><td></td><td>1.104<sub>6</sub></td><td>117*</td><td>132<sub>8</sub></td><td></td><td>3.71<sub>29</sub></td></tr><tr><td>PHTE/MXDA</td><td></td><td>1.202<sub>5</sub></td><td>176*</td><td>170<sub>10</sub></td><td></td><td>4.88<sub>18</sub></td></tr></tbody></table>

<sup>a</sup>Values with asterisks were obtained experimentally, while # indicates T<sub>α</sub>. The subscripts give the accuracy of the last decimal(s): i.e., 1.050<sub>3</sub> means 1.050 ± 0.003.

between the rapidity of the simulation and the accuracy of the results. As previously observed by Godey,<sup>54</sup> a transient region persists between the rubbery and vitreous states. This is relevant to the enlargement of the glass transition domain in a confined environment.<sup>55,56</sup> Region II is delimited by two temperatures, T<sub>1</sub> and T<sub>2</sub>, which correspond to a departure from the linear behavior of the free volume with respect to the temperature. Figure 7 shows the three regions obtained, and

<figure><img src="image_16.png" /><figcaption>Figure 7. Free volume of DGEBA+IPDA with respect to the temperature. The glass transition temperature is considered as the middle vertical line (green dashed line).</figcaption></figure>

the glass transition temperature is conventionally taken as the intercept of the linear fits of region I and III (denoted as the green dashed line).

Usually, the coefficient of volumetric thermal expansion (CVTE) is adopted to determine the linear fit of the domains. Here, this coefficient is too noisy for the thermosets studied. A piecewise linear function is thus applied to fit the three domains accordingly; the script is provided by the Python library pwlf.<sup>57</sup> This brings some automation to the fit of the linear domains, often very complicated and subjective to do. For each system, eight independent simulations are performed

and the change in free volume is measured as a function of temperature.

However, in order to compare the simulated to the experimental values, the relation of Williams-Landel-Ferry (eq 3), which relies on the principle of time-temperature equivalence, is applied by using the empirical constants C<sub>1</sub> = 17.44 and C<sub>2</sub> = 51.6 K.

$$ \log_{10} \alpha_t = \log_{10} \frac{q_{\text{exptl}}}{q_{\text{sim}}} = \frac{-C_1(T_{\text{g,shifted}} - T_{\text{g,exptl}})}{C_2 + T_{\text{g,shifted}} - T_{\text{g,exptl}}} \quad (3) $$

It was deemed acceptable to use the universal constants because of the structure of the studied resins and the time constraint of using fitted constraints for each.<sup>58</sup> In the scope of this work, several systems have been done both experimentally and by molecular simulations; this is mainly the case for the bio-based resins (DGEVA and PHTE). First, the correlation between these experimental and simulated glass transition temperatures is acceptable, with an average absolute deviation of 12%. The difference can be as minimal as 0.4% for PHTE/IPDA and as important as 59.5% for DGEVA/DETA. Second, there is a disparity between experimental results if a comparison is made between values found in this work and those from references. As a matter of fact, a mean difference of ±11 K is found, with a disparity of up to ±27 K for the case of DGEBA/EDA. Nevertheless, the results in Figure 8 show a general great correlation between simulated and experimental data. The average absolute deviation is about 7.0% without taking into account the three identified outliers and 12.2% with. Indeed, the simulated values of the glass transition temperature are well within the error range of the experiments, which confirms the interest in using molecular simulations to predict this property. The tendency is surely well-preserved even though generic WLF coefficients are applied. In cases where the prediction is less efficient, a tuning of those coefficients might be of interest. However, it will require additional time-consuming simulations of various cooling speeds for each system. The gap between experimental and simulated results regarding the DGEVA/DETA, DGEBA/
<figure><img src="image_17.png" /><figcaption>Figure 8. Comparison of simulated and experimental glass transition temperature $T_g$. Results taken from other works are given in Table 2.</figcaption></figure>

EDA, and DGEBA/MMIPA systems has to be further investigated.

**3.3. Finding Young's Modulus.** Experimentally, Young's modulus is typically obtained by traction testing. The results are widely dependent on the stress-strain applied during the experiment. Nevertheless, these types of tests are normed, and the classic speed used is on the order of mm/min. However, it is practically impossible to measure the modulus at experimental strain rates in an MD system. One way to bypass this constraint is by using an innovative method based on the fluctuation of local density⁵⁹ inside the solid epoxy materials at equilibrium. Consequently, an NVT simulation is performed during 2 ns at 300 K for each system. The structural factor S is determined from the atomistic positions and is a function of the wave vector q. The following extrapolation in eq 4 is done to obtain the bulk and shear moduli, respectively B and G. The factor 2 is added here to take into consideration the two transversal directions compared to the expression established in ref 59 for a two-dimensional system (see the Supporting Information for a demonstration of the following expressions).

$$ \lim_{q \to 0} S_{\perp}(q) = 2\rho TG \quad \lim_{q \to 0} S_{\parallel}(q) = \frac{\rho T}{B + 4/3G} \quad (4) $$

The first member of the preceding equations is obtained graphically by linearly fitting the structure factor S and acquiring the intercept, which corresponds to the limit when the wave vector q tends to zero. A visual representation of this graphical method is found in Figure 9. However, the property of interest is Young's modulus. We suppose a homogeneous isotropic linear elastic material so that its elastic properties are uniquely determined by any two moduli. Thus, having two moduli grants access to any other of the elastic moduli. The following equation is used to determine Young's modulus from the bulk and shear moduli computed previously:

$$ E = \frac{9BG}{3B + G} \quad (5) $$

Figure 10 depicts the simulated Young's modulus determined here and a few experimental points taken from a previous experiment.⁴⁸ As can be seen, the simulation fits roughly the few experimental data points, which is promising.

<figure><img src="image_18.png" /><figcaption>Figure 9. Graphical extrapolation of the structure factors.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Figure 10. Comparison of the simulated with the experimental Young's modulus. Results taken from other works are given in Table 2.</figcaption></figure>

Knowing that the simulated values exhibit a significant correlation with an average absolute deviation of 2.7%, this seems like a fair assumption. This thus raises interest in the hope that this trend will be transferable to all modeled systems. However, a wider sample is needed in terms of experimental data sets and range of Young's moduli to positively conclude.

**3.4. Influence of the Molecular Mass $M_n$.** One of the advantages of using MD simulations is to be able to straightforwardly customize the molecules used in the epoxy resins. As well as easily switching from one hardener to another, one can play with the possible oligomerization of the epoxies. The study done here enlightens the role that the molecular weight plays in the resulting thermomechanical properties of the resins. This investigation corresponds to a DGEBA+IPDA resin and includes six different oligomerizations of the DGEBA epoxy with values of n equal to 0.08, 0.32, 1, 3, 5, and 10 (see Figure 1). Values of n equal to 0.08 and 0.32 were taken from commercial references Sicomin, SR GreenPoxy 28 and Sika, Biresin CR170, respectively. Their population of n = 0, n = 1, etc. were determined by ¹H NMR. Otherwise, integer values of n were taken as a unique
<figure><img src="image_20.png" /><figcaption>**Figure 11.** Influence of the molecular weight on the thermomechanical properties of a DGEBA/IPDA system.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>**Figure 12.** Influence of the conversion rate on the thermomechanical properties of a DGEBA/IPDA system.</figcaption></figure>

An additional experimental point at $n = 5.3$ is proposed, corresponding to the Epikote 1004 commercial reference. This sample's density was experimentally measured, but a large amount of air was imprisoned inside the material, leading to an underestimation of the density.

A focus is made on the influence of the molecular mass on density, glass transition temperature, and Young's modulus (see Figure 11). Previous studies on DGEBA and Jeffamines show the link between the length of the cross-linker and the resulting properties.39,40 Similarly, this relationship is also displayed here: an increase in molecular mass causes an increase of the density and a decrease of the glass transition temperature and Young's modulus. This outcome can be explained first by the flexibility that the increase of the distance between cross-linking points brings, having a higher degree of oligomerization of the epoxies. By having this oligomerization, it is then easier to find a more relaxed structure, thus bringing the molecules closer together and increasing the density. The length of a cross-link is often one of the main focuses when it comes to tuning thermomechanical properties of materials. One could want to limit the brittleness of a material by lengthening a cross-link to bring some flexibility to a material, but it needs to be taken into account that this change is typically detrimental to the thermal stability. This study enlightens the dual effect of the cross-link length and how a middle ground needs to be reached when designing a new material, and how molecular simulations can help to that effect.

**3.5. Influence of the Conversion Rate.** In the same spirit, a study on the influence of the conversion rate was done to enlighten another facet of the complexity of the structure–property relationship (Figure 12). In the scope of this work, a DGEBA/IPDA resin was studied with various curing rates of 30%, 50%, 70%, and 90%. When it comes to conversion rate, this link has been shown previously on Jeffamine systems.39 First, a slight decrease in density was observed here while it was not the case for the Jeffamine systems. This could be explained by the fact that an increase in cross-links could further constrain the system. The limited opportunity for monomers to be optimally stacked could lead to an increase in free volume and thus result in a decrease in density. Second, the trend regarding the glass transition temperature and Young's modulus is coherent with that of the Jeffamine systems. The link between the structure and the resulting properties is more readily accessible than for the density observations. Since a low conversion rate accounts for fewer cross-links, it is coherent to have a higher elastic response for this material. Similarly, less thermal energy is needed to go from a rigid amorphous thermoset to a more flexible one when fewer covalent bonds are present.

## 4. CONCLUSION

This work enlightens a broader study of the modeling of high-performance materials such as epoxy thermosets in terms of the chemical structure and thermomechanical properties. Novel resins were investigated, including bio-based ones,
both experimentally and by molecular dynamics. The implementation of a “reactive” cross-linking algorithm allows for the creation of a highly cross-linked resin, of up to 99%. The use of the CHARMM force field, usually reserved for biomolecules, was proven to be suitable for these thermosetting resins. Indeed, this work exhibits a good correlation between experimental and simulated results for a variety of resins with backbones fluctuating from aliphatic to aromatic. Thermomechanical properties in the name of density, Young's modulus, and glass transition temperature were measured. In order to acquire these properties by simulation, innovative methodologies and post-treatment of the data were performed and compared with the experimental data to fully implement a robust and generalized approach for each property studied. The most compelling comparison is found in the glass transition temperature, as it is the property with the most data points both experimentally and by molecular simulation. Also, the novelty of this work resides in the large range of glass transition temperatures covered by the resins studied. Overall, because of the consistency with the experimental data, with average absolute deviations between experiments and molecular simulation comprised between 2% and 12%, it is possible to conclude that all of the methods presented are transferable and successful. Thus, it is possible to efficiently determine the density, the glass transition temperature, and Young's modulus for epoxy thermosets by atomistic methods. From an industrial point of view, the development of these methods is the starting point of a more ambitious goal. The definition of these first properties allows a quick characterization of a system and an evaluation if it is interesting for an experimental study according to the targeted specifications. As well as being a predicting tool, this investigation demonstrates the interest in using molecular simulations to extract important structural analyses, often difficult to obtain experimentally. Atomistic simulations will thus have a crucial role in the structure—property understanding of materials by visualizing the system and acquiring key properties. Building from this, other types of resins such as polyurethanes are currently being considered and studied with this approach.

## ASSOCIATED CONTENT

### Supporting Information

The Supporting Information is available free of charge at
https://pubs.acs.org/doi/10.1021/acsomega.2c03071.

Demonstration of the operational expressions of $S_⊥$ and $S_‖$ in a three-dimensional system (PDF)

## AUTHOR INFORMATION

### Corresponding Authors

Mathilde Orselly – *Specific Polymers, 34160 Castries, France; Université Clermont Auvergne, Clermont Auvergne INP, CNRS, Institut de Chimie de Clermont-Ferrand, F-63000 Clermont-Ferrand, France*; Email: mathilde.orselly@specificpolymers.fr

Patrice Malfreyt – *Université Clermont Auvergne, Clermont Auvergne INP, CNRS, Institut de Chimie de Clermont-Ferrand, F-63000 Clermont-Ferrand, France*; orcid.org/0000-0002-3710-5418; Email: patrice.malfreyt@uca.fr

## Authors

Julien Devemy – *Université Clermont Auvergne, Clermont Auvergne INP, CNRS, Institut de Chimie de Clermont-Ferrand, F-63000 Clermont-Ferrand, France*

Agathe Bouvet-Marchand – *Specific Polymers, 34160 Castries, France*; orcid.org/0000-0003-0371-1126

Alain Dequidt – *Université Clermont Auvergne, Clermont Auvergne INP, CNRS, Institut de Chimie de Clermont-Ferrand, F-63000 Clermont-Ferrand, France*; orcid.org/0000-0003-1206-1911

Cédric Loubat – *Specific Polymers, 34160 Castries, France*

Complete contact information is available at:
https://pubs.acs.org/10.1021/acsomega.2c03071

## Notes

The authors declare no competing financial interest.

## ACKNOWLEDGMENTS

M.O., A.B.-M. and C.L. thank all the members of SimatLab for stimulating discussions about this work. SimatLab is a joint public-private laboratory dedicated to the multiscale modeling of polymer materials. This laboratory is supported by Michelin, Clermont Auvergne University (UCA), CHU of Clermont-Ferrand, and the CNRS. We are also grateful to the Mésocentre Clermont Auvergne University for providing computing and storage resources. We would like to thank Professor L. Berthier for fruitful discussions on the derivation of Equation 4.

## REFERENCES

(1) Pascault, J.-P.; Williams, R. J. J. *Epoxy Polymers*; Wiley: 2010; Chapter 1, pp 1–12.

(2) Li, C.; Strachan, A. *Molecular Scale Simulations on Thermoset Polymers: A Review*. *J. Polym. Sci. B Polym. Phys.* **2015**, *53*, 103–122.

(3) Komarov, P. V.; Yu-Tsung, C.; Shih-Ming, C.; Khalatur, P. G.; Reineker, P. Highly Cross-Linked Epoxy Resins: An Atomistic Molecular Dynamics Simulation Combined with a Mapping/Reverse Mapping Procedure. *Macromolecules* **2007**, *40*, 8104–8113.

(4) Gavriloğlu, A. A.; Komarov, P. V.; Khalatur, P. G. Thermal Properties and Topology of Epoxy Networks: A Multiscale Simulation Methodology. *Macromolecules* **2015**, *48*, 206–212.

(5) Guseva, D. V.; Rudyak, V. Y.; Komarov, P. V.; Sulimov, A. V.; Bulgakov, B. A.; Chertovich, A. V. Crosslinking Mechanisms, Structure and Glass Transition in Phthalonitrile Resins: Insight from Computer Multiscale Simulations and Experiments. *J. Polym. Sci. B Polym. Phys.* **2018**, *56*, 362–374.

(6) Prasad, A.; Grover, T.; Basu, S. Coarse – Grained Molecular Dynamics Simulation of Cross-linking of DGEBA Epoxy Resin and Estimation of the Adhesive Strength. *Int. J. Eng. Sci. Technol.* **2010**, *2*, 1.

(7) Gavrielides, A.; Duguet, T.; Aufray, M.; Lacaze-Dufaure, C. Model of the DGEBA-EDA Epoxy Polymer: Experiments and Simulation Using Classical Molecular Dynamics. *Int. J. Polym. Sci. 2019*, *2019*, 1.

(8) Radue, M. S.; Varshney, V.; Baur, J. W.; Roy, A. K.; Odegard, G. M. Molecular Modeling of Cross-Linked Polymers with Complex Cure Pathways: A Case Study of Bismaleimide Resins. *Macromolecules* **2018**, *51*, 1830–1840.

(9) Bandyopadhyay, A.; Valavala, P. K.; Clancy, T. C.; Wise, K. E.; Odegard, G. M. Molecular Modeling of Crosslinked Epoxy Polymers: The Effect of Crosslink Density on Thermomechanical Properties. *Polymer* **2011**, *52*, 2445–2452.

(10) Fan, H. B.; Yuen, M. M. Material Properties of the Cross-linked Epoxy Resin Compound Predicted by Molecular Dynamics Simulation. *Polymer* **2007**, *48*, 2174–2178.
(11) Nouri, N.; Ziaei-Rad, S. A Molecular Dynamics Investigation on Mechanical Properties of Cross-Linked Polymer Networks. *Macromolecules* **2011**, *44*, 5481–5489.

(12) Li, C.; Strachan, A. Molecular simulations of crosslinking process of thermosetting polymers. *Polymer* **2010**, *51*, 6058–6070.

(13) Shokuhfar, A.; Arab, B. The Effect of Cross-linking Density on the Mechanical Properties and Structure of the Epoxy Polymers: Molecular Dynamics Simulation. *J. Mol. Model.* **2013**, *19*, 3719–3731.

(14) Gissinger, J. R.; Jensen, B. D.; Wise, K. E. Modeling chemical reactions in classical molecular dynamics simulations. *Polymer* **2017**, *128*, 211–217.

(15) Gissinger, J. R.; Jensen, B. D.; Wise, K. E. REACTER: AHeuristic Method for Reactive Molecular Dynamics. *Macromolecules* **2020**, *53*, 9953–9961.

(16) Jang, C.; Sirk, T. W.; Andzelm, J. W.; Abrams, C. F. Comparison of Crosslinking Algorithms in Molecular Dynamics Simulation of Thermosetting Polymers. *Macromol. Theory Simul.* **2015**, *24*, 260–270.

(17) Tack, J. L.; Ford, D. M. Thermodynamic and Mechanical Properties of Epoxy Resin DGEBF Crosslinked with DETDA by Molecular Dynamics. *J. Mol. Graph. Model.* **2008**, *26*, 1269–1275.

(18) Varshney, V.; Patnaik, S. S.; Roy, A. K.; Farmer, B. L. A Molecular Dynamics Study of Epoxy-Based Networks: Cross-Linking Procedure and Prediction of Molecular and Material Properties. *Macromolecules* **2008**, *41*, 6837–6842.

(19) Yu, S.; Yang, S.; Cho, M. Multiscale Modeling of Cross-linked Epoxy Nanocomposites to Characterize the Effect of Particle Size on Thermal Conductivity. *J. Appl. Phys.* **2011**, *110*, 124302.

(20) Choi, J.; Yu, S.; Yang, S.; Cho, M. The Glass Transition and Thermoelastic Behavior of Epoxy-based ?anocomposites: A Molecular Dynamics Study. *Polymer* **2011**, *52*, 5197–5203.

(21) Yu, S.; Yang, S.; Cho, M. Multi-scale Modeling of Cross-linked Epoxy Nanocomposites. *Polymer* **2009**, *50*, 945–952.

(22) Schichtel, J. J.; Chattopadhyay, A. Modeling Thermoset Polymers Using an Improved Molecular Dynamics Crosslinking Methodology. *Comput. Mater. Sci.* **2020**, *174*, 109469.

(23) Wu, C.; Xu, W. Atomistic Molecular Modelling of Crosslinked Epoxy Resin. *Polymer* **2006**, *47*, 6004–6009.

(24) Wu, C.; Xu, W. Atomistic Molecular Simulations of Structure and Dynamics of Crosslinked Epoxy Resin. *Polymer* **2007**, *48*, 5802–5812.

(25) Wang, Z.; Lv, Q.; Chen, S.; Li, C.; Sun, S.; Hu, S. Glass Transition Investigations on Highly Crosslinked Epoxy Resins by Molecular Dynamics Simulations. *Mol. Simul.* **2015**, *41*, 1515–1527.

(26) Lin, P.-H.; Khare, R. Molecular Simulation of Cross-Linked Epoxy and Epoxy-POSS Nanocomposite. *Macromolecules* **2009**, *42*, 4319–4327.

(27) Jeyranpour, F.; Alahyarizadeh, G.; Arab, B. Comparative Investigation of Thermal and Mechanical Properties of Cross-linked Epoxy Polymers with Different Curing Agents by Molecular Dynamics Simulation. *J. Mol. Graph. Model.* **2015**, *62*, 157–164.

(28) Clancy, T.; Frankland, S.; Hinkley, J.; Gates, T. Molecular Modeling for Calculation of Mechanical Properties of Epoxies with Moisture Ingress. *Polymer* **2009**, *50*, 2736–2742.

(29) Demir, B.; Walsh, T. R. A Robust and Reproducible Procedure for Cross-linking Thermoset Polymers using Molecular Simulation. *Soft Matter* **2016**, *12*, 2453–2464.

(30) Liu, H.; Li, M.; Lu, Z.-Y.; Zhang, Z.-G.; Sun, C.-C.; Cui, T. Multiscale Simulation Study on the Curing Reaction and the Network Structure in a Typical Epoxy System. *Macromolecules* **2011**, *44*, 8650–8660.

(31) Chang, S.-H.; Kim, H.-S. Investigation of Hygroscopic Properties in Electronic Packages Using Molecular Dynamics Simulation. *Polymer* **2011**, *52*, 3437–3442.

(32) Soni, N. J.; Lin, P.-H.; Khare, R. Effect of Cross-linker Length on the Thermal and Volumetric Properties of Cross-linked Epoxy Networks: A Molecular Simulation Study. *Polymer* **2012**, *53*, 1015–1019.

(33) Sirk, T. W.; Karim, M.; Khare, K. S.; Lenhart, J. L.; Andzelm, J. W.; Khare, R. Bi-modal Polymer Networks: Composition-Dependent Trends in Thermal, Volumetric and Structural Properties from Molecular Dynamics Simulation. *Polymer* **2015**, *58*, 199–208.

(34) Arab, B.; Shokuhfar, A.; Ebrahimi-Nejad, S. Glass Transition Temperature of Cross-linked Epoxy Polymers: a Molecular Dynamics Study. *Proceedings of the International Conference Nanomaterials: Applications and Properties*, 2012; pp 1–4.

(35) Arab, B.; Shokuhfar, A. Molecular Dynamics Simulation of Cross-Linked Epoxy Polymers: the Effect of Force Field on the Estimation of Properties. *J. Nano-Electron. Phys.* **2013**, *5*, 01013.

(36) Langeloth, M.; Sugii, T.; Böhm, M. C.; Müller-Plathe, F. The Glass Transition in Cured Epoxy Thermosets: A Comparative Molecular Dynamics Study in Coarse-Grained and Atomistic Resolution. *J. Chem. Phys.* **2015**, *143*, 243158.

(37) Masoumi, S.; Arab, B.; Valipour, H. A study of Thermo-mechanical Properties of the Cross-linked Epoxy: An Atomistic Simulation. *Polymer* **2015**, *70*, 351–360.

(38) Fan, J.; Anastassiou, A.; Macosko, C. W.; Tadmor, E. B. Molecular Dynamics Predictions of Thermomechanical Properties of an Epoxy Thermosetting Polymer. *Polymer* **2020**, *196*, 122477.

(39) Chowdhury, S. C.; Elder, R. M.; Sirk, T. W.; Gillespie, J. W. Epoxy Resin Thermo-mechanics and Failure Modes: Effects of Cure and Cross-linker Length. *Compos. B. Eng.* **2020**, *186*, 107814.

(40) Patil, P. N.; Rath, S. K.; Sharma, S. K.; Sudarshan, K.; Maheshwari, P.; Patri, M.; Praveen, S.; Khandelwal, P.; Pujari, P. K. Free Volumes and Structural Relaxations in Diglycidyl Ether of Bisphenol-A Based Epoxy-Polyether Amine Networks. *Soft Matter* **2013**, *9*, 3589.

(41) Vanommeslaeghe, K.; Hatcher, E.; Acharya, C.; Kundu, S.; Zhong, S.; Shim, J.; Darian, E.; Guvench, O.; Lopes, P.; Vorobyov, I.; Mackerell, A. D. CHARMM General Force Field: A Force Field for Drug-Like Molecules Compatible with the CHARMM All-Atom Additive Biological Force Fields. *J. Comput. Chem.* **2009**, *31*, 671–690.

(42) Vanommeslaeghe, K.; MacKerell, A. D. Automation of the CHARMM General Force Field (CGenFF) I: Bond Perception and Atom Typing. *J. Chem. Inf. Model.* **2012**, *52*, 3144–3154.

(43) Vanommeslaeghe, K.; Raman, E. P.; MacKerell, A. D. Automation of the CHARMM General Force Field (CGenFF) II: Assignment of Bonded Parameters and Partial Atomic Charges. *J. Chem. Inf. Model.* **2012**, *52*, 3155–3168.

(44) Martínez, L.; Andrade, R.; Birgin, E. G.; Martínez, J. M. PACKMOL: A Package for Building Initial Configurations for Molecular Dynamics Simulations. *J. Comput. Chem.* **2009**, *30*, 2157–2164.

(45) Nosé, S. A Molecular Dynamics Method for Simulations in the Canonical Ensemble. *Mol. Phys.* **1984**, *52*, 255–268.

(46) Plimpton, S. Fast Parallel Algorithms for Short-range Molecular Dynamics. *J. Comput. Chem.* **1995**, *117*, 1–19.

(47) Ernault, E.; Richaud, E.; Fayolle, B. Thermal-Oxidation of Epoxy/Amine Followed by Glass Transition Temperature changes. *Polym. Degrad. Stab.* **2017**, *138*, 82–90.

(48) Garcia, F. G.; Soares, B. G.; Pita, V. J. R. R.; Sanchez, R.; Rieumont, J. Mechanical Properties of Epoxy Networks Based on DGEBA and Aliphatic Amines. *J. Appl. Polym. Sci.* **2007**, *106*, 2047–2055.

(49) Aufray, M.; Roche, A. A. Residual Stresses and Practical Adhesion: Effect of Organo-Metallic Complex Formation and Crystallization. *J. Adhes. Sci. Technol.* **2006**, *20*, 1889–1903.

(50) Jilani, W.; Mzabi, N.; Fourati, N.; Zerrouki, C.; Gallot-Lavallée, O.; Zerrouki, R.; Guermazi, H. Effects of Curing Agent on Conductivity, Structural and Dielectric Properties of an Epoxy Polymer. *Polymer* **2015**, *79*, 73–81.

(51) Prolongo, S.; Mikes, F.; Cabanelas, J.; Paz-Abuin, S.; Baselga, J. Diffusion Control on the Cure Kinetics of DGEBA with Ethylenediamines. *J. Mater. Process. Technol.* **2003**, *143–144*, 546–550.
(52) Fache, M.; Auvergne, R.; Boutevin, B.; Caillol, S. New Vanillin-Derived Diepoxy Monomers For the Synthesis of Biobased Thermosets. *Eur. Polym. J.* **2015**, *67*, 527–538.

(53) Pascault, J.-P.; Sautereau, H.; Verdu, J.; Williams, R. J. J. *Thermosetting Polymers*; CRC Press: 2002.

(54) Godey, F.; Bensaid, M.; Soldera, A. Extent of the Glass Transition in Polymers Envisioned by Computation of Mechanical Properties. *Polymer* **2019**, *164*, 33–38.

(55) Vignaud, G.; S. Chebil, M.; Bal, J. K.; Delorme, N.; Beuvier, T.; Grohens, Y.; Gibaud, A. Densification and Depression in Glass Transition Temperature in Polystyrene Thin Films. *Langmuir* **2014**, *30*, 11599–11608.

(56) Lan, T.; Torkelson, J. M. Methacrylate-Based Polymer Films Useful in Lithographic Applications Exhibit Different Glass Transition Temperature-Confinement Effects at High and Low Molecular Weight. *Polymer* **2014**, *55*, 1249–1258.

(57) Jekel, C. F.; Venter, G. pwlf: A Python Library for Fitting 1D Continuous Piecewise Linear Functions, 2019.

(58) Dudowicz, J.; Douglas, J. F.; Freed, K. F. The Meaning of the “Universal” WLF Parameters of Glass-Forming Polymer Liquids. *J. Chem. Phys.* **2015**, *142*, 014905.

(59) Klix, C. L.; Ebert, F.; Weysser, F.; Fuchs, M.; Maret, G.; Keim, P. Glass Elasticity From Particle Trajectories. *Phys. Rev. Lett.* **2012**, *109*, 178301.