Model of the DGEBA-EDA Epoxy Polymer: Experiments and Simulation Using Classical Molecular Dynamics

Polyepoxy samples are synthesized from diglycidylether of bisphenol A (DGEBA) and ethylene diamine (EDA) monomers at a stoichiometric ratio of 2 DGEBA:1 EDA in model conditions in order to promote a high degree of polymerization and a low density of defects and to try to approach the ideal models obtained by simulation. A slow polymerization ( > 24 h at ambient temperature) and a postcuring achieved in an inert atmosphere lead to a conversion degree of 92 ± 2% and a midpoint glass transition temperature of 391 ± 1K . In parallel, a model is created with a multistep cross-linking procedure. In this work, all-atom molecular dynamics (MD) simulations are performed with LAMMPS and the GAFF 1.8 force ﬁ eld. In the initial liquid mixture of reactants (600 molecules), proper mixing is demonstrated by the calculation of the partial radial distribution functions (RDF), which show a minimum intermolecular distance of 2.8Å and similar distributions for EDA-EDA, DGEBA-DGEBA, and DGEBA-EDA molecules in the simulation boxes. Then, in alternation with MD equilibrations, cross-linking is performed on frozen con ﬁ gurations by creating covalent bonds between reactive pairs within a reaction radius of 3Å. The resulting boxes show conversion rates of 90-93% and densities close to the experimental value. Finally, a cooling ramp from 700 K to 25 K is applied in order to monitor the speci ﬁ c volume and the coe ﬃ cient of volumetric thermal expansion (CVTE) of the polymer and to derive the glass transition temperature. Experimental thermomechanical analyses (TMA) compares well with simulations for both the speci ﬁ c volume and the CVTE evolutions with temperature. Whereas the uncertainty remains high with the ﬁ tting procedure used, we calculate a glass transition temperature of 390 ± 8 K which compares very well with the experimental values ( 391 ± 1K from DSC and 380 K from TMA).


Introduction
One of the great challenges in computational chemistry is to derive the physical, mechanical, and thermodynamic properties of solids and especially of polymers. In recent years, polymers have been implemented in a variety of applications, from medicine to space industries [1][2][3][4][5]. Therefore, it is of great importance to investigate their network structures and properties in order to engineer their bulk and surface properties efficiently.
Thermosetting polymers of the polyepoxy family are widely used in applications in the areas of aeronautics and space, luxury, and sports. They exhibit superior thermal and mechanical properties when compared to other polymers, and thus are used as the matrix in composite materials (e.g., carbon fiber/epoxy). Unfortunately, such composites lack a variety of properties desirable for many components, such as electrical and thermal conductivities. Surface functionalization through metallization is then necessary, but it still represents a challenge [6,7]. A better knowledge of polyepoxy surfaces would be a first step towards the understanding of surface phenomena occurring during metallization.
In the present work, we use atomistic MD and a static cross-linking code to create a model polyepoxy from the epoxy resin diglycidyl ether of bisphenol-A (DGEBA) and ethylenediamine (EDA).
We have been studying this model system in the last few years because the monomers are small and suitable for model experiments and ab initio calculations [30,31], but yet close to real materials used in the industry. To generate atomistic model configurations, we have developed a cross-linking procedure based on the works of Jang et al. [27] and Sirk et al. [32]. Molecular dynamics simulations at 700 K are used in order to mix efficiently and relax the strain in the simulation boxes and thus achieve the maximum percentage of polymerization starting from a stoichiometric mixture of DGEBA and EDA (2 : 1). The thermal properties of the model polyepoxy are characterized by the determination of the density, the coefficient of volumetric thermal expansion, and the glass transition temperature (T g ) of the polymer, which are compared to experimental data.

Experimental Details.
The chemical structures of the epoxy monomer (DGEBA) and the hardener (the primary diamine EDA) are shown in Figure 1 (reactants 1 and 2), along with a few polyaddition steps of the polymerization reaction. DGEBA is a bifunctional reactant with two epoxide groups, while EDA has four reactive sites with two primary amine groups. Thus, the stoichiometry ratio of DGE-BA : EDA is 2 : 1. The polymerization consists in the reaction of the H atoms of the amine groups of EDA with the O atoms of the epoxide groups of DGEBA and the bonding of the N atoms of the amine group with the C atoms of the epoxide groups. We assume that the homopolymerization is negligible in the present conditions of temperature and in the absence of a catalyst.
We use the stoichiometric ratio of DGEBA (D.E.R. 332, Dow Chemical Co., n = 0 03) and EDA (Analytical Grade, purity > 99.5%, Sigma-Aldrich). Since the mass of DGEBA (m DGEBA ) is fixed to 5 g, the mass of EDA (m EDA ) is determined as follows: where M DGEBA is the molar mass (348.52 g/mol) of this DGEBA and f DGEBA is its functionality (2), and M EDA is the molar mass (60.10 g/mol) of EDA and f EDA is its functionality (4). We assume that no etherification occurs. The DGEBA resin has a density of 1.170 g/cm 3 , and the onset of its glass transition temperature determined by differential scanning calorimetry equals 231 K. EDA is a liquid amine with a density of 0.897 g/cm 3 . The mixture is then mechanically stirred in an Ar glove box for 7 min before it is poured as a thin droplet on Si coupons. Polymerization is then allowed for 48 h at ambient temperature, followed by a postcuring of 2 h at 413 K, either in the Ar glovebox, in ambient atmosphere, or in a primary vacuum (with a short air exposure in this latter case).
Fourier transform infrared spectroscopy (FTIR, Frontier, PerkinElmer equipped with a NIR TGS detector) is achieved in transmission in the 4000-8000 cm -1 range. 16 scans are collected for each analysis with a resolution of 4 cm -1 . Background spectra are performed on the two empty microscope glass slides used for maintaining the polyepoxy samples. The decrease in the IR normalized band intensity of epoxy groups is directly proportional to the increase of conversion degree. We quantify the characteristic epoxy band (combination band of the -CH 2 of the epoxy group) at 4530 cm -1 . The reference band (i.e., the band conserving the same area during the chemical reaction of polymerization) is the combination band of C=C with aromatic -CH at 4623 cm -1 [33]. Peak areas are then used for calculating the conversion degree (X e,NIR ) of epoxy groups: The terms A epoxy and A reference refer to the areas of the peaks of the epoxy and reference groups, respectively (i.e., A epoxy = A 4530 and A ref = A 4623 ).
Differential scanning calorimetry (DSC) is used for the determination of the glass transition temperature (T g ). We use a DSC 204 Phoenix Series (NETZSCH) coupled with a TASC 414/4 controller. The apparatus is calibrated against melting temperatures of In, Hg, Sn, Bi, and Zn, applying a +10 deg/min temperature ramp. Samples are placed in aluminum capsules. Mass is measured with an accuracy of ±0.1 mg. We report and compare the midpoint T g temperatures.

2
International Journal of Polymer Science In addition, thermomechanical analyses (TMA) are performed on a cubic 3 7 × 3 7 × 3 7 mm 3 sample with a NETZSCH TMA 402 F3 apparatus. 3 cycles are executed consecutively from 134 K to 485 K with a heating ramp of +5 deg/min. During the first cycle, the sample is heated up and maintained at 485 K for 1 h before cooling. This is intended to erase the thermal history of the sample. The other two cycles include a 15 min plateau only at maximum temperature.

Computational Details.
The generalized AMBER force field (GAFF 1.8) [34] is used to describe the intra-and intermolecular interactions. Any missing force field parameters are determined using the AmberTools16 package [35]. The energy minimizations and molecular dynamics calculations are carried out using the LAMMPS simulation package [36]. The atomic charges on all atoms were calculated with the RESP method [37], applied to monomers, dimers (1 DGEBA bonded to 1 EDA), and trimers (2 DGEBA bonded to 1 EDA) in order to get charges for all atom types present during the cross-linking procedure and in the final polyepoxy. The particle particle particle mesh (PPPM) algorithm [38] is used to calculate the coulombic interactions with a real-space cutoff of 9 Å. A cutoff distance of 9 Å (with the buffer of 2 Å) is also used for the van der Waals interactions.
The velocity-Verlet integrator [39] is used to integrate the equations of motion using a time step of 1 fs. The temperature and pressure are controlled using the Nosé-Hoover thermostat [40] and barostat [35].
The first step of our methodology is to generate a "mixture of reactants" containing the stoichiometric amounts of DGEBA and EDA. In order to approximate a correct representation of a bulk polyepoxy, we build a cubic box which contains 400 DGEBA and 200 EDA molecules at a density fixed at 0.7 g/cm 3 using PACKMOL [41] and the amorphous builder module included in the Materials Studio package [42].
The process of obtaining a polymer from the initial monomer mixture involves the cross-linking of monomers. To that aim, we developed a homemade code (FORTRAN 90) based on the works of Jang et al. [27] and Sirk et al. [32]. This multistep cross-linking process is presented in Figure 2. MD simulations are performed at 700 K between two cross-linking steps in order to relax the structures and decrease local strain. In complement, two runs are performed at the experimental curing temperature (413 K) and at high temperature (900 K).
The multistep cross-linking process sketched in Figure 2 consists in the following steps: (1) Packing of the reactants (stoichiometric mixture: 400 DGEBA and 200 EDA molecules) in a box of predefined density (0.7 g/cm 3 ) (2) NVE MD simulation (constant number of atoms, volume, and energy) for 1 ns, followed by NPT simulation (constant number of atoms, pressure, and temperature) for 1 ns, at a pressure P = 1 atm and at a temperature T = 700 K to equilibrate the initial monomer mixture 3 International Journal of Polymer Science (3+4) Pairs of unreacted atoms are identified within a search reaction cutoff distance of 3 Å. To cross-link these atoms, the identities (IDs) of modified atoms are changed. The topology file is updated. New angles, dihedrals, and RESP charges that result from the creation of the new cross-linked bonds are assigned (5) MD simulations (1 ns) in NVE, NVT, and NPT (700 K, 1 atm) ensembles are performed; the temperature is fixed at 700 K, which is far above the glass transition temperature T g of the polymer determined experimentally (391 K), thus allowing some mobility of the chains to obtain a relaxed structure The sequence of (3 + 4 + 5) steps is repeated until no change in the conversion degree is observed between two sequences. The maximal conversion degree is then obtained.
(6) An NPT MD simulation at 700 K for 2 ns is performed on the final reticulated simulation box, followed by stepwise cooling to 300 K to get relaxed structures and density of the polyepoxy model at 300 K. More details on the cooling process are given below in Section 3.2.4 where results on the volume-temperature behavior of the polymer are presented It is noted that the cross-linking algorithm outlined above is not intended to mimic the actual polymerization reaction kinetics, as conducted in the laboratory. Experimentally, curing is carried out at temperatures of 140°C (413 K) over much longer times than can be accessed by atomistic MD simulation. Our choice of 700 K, a temperature actually exceeding the thermal decomposition temperature of ordinary cured epoxy resins, is motivated by the need to bring about a vigorous configurational rearrangement within the time spans of tens of nanoseconds that can be simulated atomistically with MD and thus push unreacted reactive groups to come together, leading to a degree of conversion comparable to that obtained experimentally. Our use of such a high temperature also accelerates the process of relaxation of the structure around a newly created N-C chemical bond. The cutoff distance of 3 Å used in the search of pairs of atoms to be bonded by chemical reaction is another parameter of the computational cross-linking procedure that does not have a direct chemical counterpart. It will thus be an interesting question whether the structural, volumetric, and thermal properties of the resulting cross-linked polymer match those obtained experimentally.

Experimental
Results. We perform experimental tests in different postcuring conditions with the aim to maximize the conversion degree. Each test is repeated five times. The results are summarized in Table 1, where the conversion degrees are derived from FTIR analyses.
Under noninert conditions, the conversion degree is low at70 ± 6%. It is improved under vacuum with a conversion degree of 84 ± 4%. However, the postcuring in the Ar glovebox ensures the highest conversion degree of 92 ± 2% and the repeatability is excellent. These results are in agreement with previous results with the DGEBA-EDA system [30].
We determine the glass transition temperature by DSC. The midpoint temperature is 391 ± 1 K, in good agreement with values reported elsewhere [43,44].   Figure 3 and they help in determining if the mixture is homogeneous. We first observe that there are no intermolecular distances below 2.8 Å. Then, except for differences between RDF amplitudes, distributions are similar and oscillate around g r = 1 (dotted line), corresponding to a homogeneous distribution. This point is important to ensure that epoxide groups are in the vicinity of amine groups to promote cross-linking and mimic the effect of the experimental mechanical stirring. We also note that the highest amplitude of the partial RDFs is when EDA molecules are present, likely because of the rather high temperature for the amine.

Polymerization
Procedure. Five boxes are generated and equilibrated as mentioned above. All boxes are then polymerized with the multistep procedure achieved at 700 K as described in Figure 2, and the conversion degree (number of created -OH over the total number of epoxy groups: N OH /2N DGEBA ) is calculated at each step and plotted in Figure 4. We observe that final conversion degrees are all above 90% after 11 steps. In the most promising cases, we obtain a conversion degree of 93%. The mobility of the remaining molecules or the steric availability of the remaining reactive sites is too low at these polymerization rates to promote further cross-linking. We also tried to maximize the polymerization rate with an additional box equilibrated in NVT and NPT at 900 K, but without more success (92%). This tends to show that the 93% limit will be difficult to surpass with these simulation parameters. As a complement to our comment in the last paragraph of Section 2.2, we show the results obtained with the multistep procedure achieved at 413 K; i.e., the temperature that is used experimentally. A plateau is reached at about 70% conversion (72% after 16 steps, not shown). This confirms that higher temperatures must be chosen if one want to achieve high conversion degrees for such short simulations times. Finally, the first conversion degrees for given steps are different from one box to another; for instance at step 2, the conversion degree of Box#1 is 70%, whereas for Box#5 it is 48%. It is noted that there was no noticeable difference between their RDFs before polymerization. Box#1 was generated in PACK-MOL, whereas Box#2-5 were generated with the randomization algorithm of Materials Studio. Since randomization   5 International Journal of Polymer Science relies on Monte Carlo algorithms, each box distribution is different and may take different trajectories during MD equilibration. Anyway, it seems that whatever the initial box configuration, the final conversion rates converge to 90% and a little more.
Box#4 is selected in its final configuration (conversion degree = 93%) in order to analyze its bonded and nonbonded RDF functions. Results are shown in Figure 5 along with a 3D representation of the box in the inset. The minimum nonbonded distances found in the polymer correspond to a peak ranging from 2 to 3 Å, i.e., the packing is dense as expected for cross-linked epoxy networks, where the interchain distance is constrained by reticulation. The nonbonded RDF rapidly reaches a value of 1 (between r = 3 5 to 5.5 Å) that demonstrates the homogeneity of the distribution of the disordered structure.
On the bonded RDF in Figure 5, we identify 5 local maxima that correspond to first-neighbor distances. Labels 1 to 5 correspond to (1)

Density.
The polymer density, calculated at 300 K and 1 atm, is also in very good concordance with the literature. We calculate a density of 1 118 ± 0 003 g/cm 3 , whereas the value obtained through empirical calculations is 1.142 g/cm 3 . This empirical density is calculated through the relationship ρ = 350 + 120M (kg/m 3 ), where M is the molecular weight of the repeat unit [44,45]. The value of M is representative of the polyepoxy polymer based on the 2DGEBA : 1EDA stoichiometry.

Glass Transition Temperature and Coefficient of
Volumetric Thermal Expansion. Experimentally, the TMA analysis consists in the monitoring of the linear expansion of the polyepoxy sample during a thermal cycle from 134 K to 485 K. We then extrapolate the measured Δl/l into a volumic Δv/v, considering an isotropic expansion. Then, we plot the specific volume v (the inverse of the calculated density) as a function of temperature. Results are shown with (black) full circles in Figure 6. For calculations, we follow a protocol based on the work of Sirk et al. [32]. First, the procedure consists in heating the polymer at 700 K for 2 ns and cooling it at a rate of 25 K/ns down to 25 K (red, empty circles in Figure 6). Two different regions are observed below and above 400 K, roughly, corresponding to the vitreous and rubbery states.
The determination of T g is explained and involves the use of the coefficient of volumetric thermal expansion (CVTE).
An elegant procedure to analyze the plots of Figure 6 consists in distinguishing the glassy and rubbery temperature domains. Actually, although there is a change in physical properties during the glass transition, the structure evolves through an undefined range of temperatures. Temperature regions must then be identified to fit the linear regimes properly and to estimate T g at the intersection. For that reason, we plot the CVTE versus temperature in Figure 7 for both the experimental and the calculation procedures.
CVTEs are calculated as follows: The term ∂v/∂T P is calculated using a finite-difference method in Origin software.
Since glass transition is a second order transition, the CVTE should show a smooth step separating two constant regions. Figure 7(a) shows the results for 3 different simulation boxes (empty circles) and for one TMA cycle (full black circles). Figure 7(b) makes a focus on the 3 TMA cycles in order to extract quantitative data more easily.
By simulation (Figure 7(a)), the glassy state region is easily identified up to T = 275 K and the dispersion between the 3 v t datasets is small. There, the average CVTE varies from 2 0 × 10 −4 K −1 at 100 K to 2 5 × 10 −4 K −1 at 275 K. At  International Journal of Polymer Science higher temperatures, the CVTE reaches 8 × 10 −4 K −1 and more, but the dispersion is high and no region can easily be identified. Nevertheless, we tentatively place a second limit at T = 450 K. By doing so, we assume that all points above 575 K are not valid because they exceed known limits for the CVTE of thermosets in the rubbery region, i.e., 4 × 10 −4 − 8 × 10 −4 K −1 [45]. An additional trial has been made to double-check the accuracy of the MD calculations at high temperatures. We applied a cooling ramp of 5 K/ns instead of 25 K/ns from 700 to 600 K, in order to equilibrate the box more smoothly (not shown). This did not improve the results, neither for the linear fitting of v T nor for the determination of a rubbery state plateau below 8 × 10 −4 K −1 . Consequences are that we can fit the simulated plot in Figure 6 with a single linear regression calculated from 25 K to 275 K. Above, there are multiple possibilities for the range of the linear regression, resulting in a high uncertainty of ±8 degrees in the determination of T g . With the linear fittings proposed, we determine a simulated T g of 390 ± 8 K. Experimentally (Figure 7(b)), the CVTE evolves linearly in the glassy region from 0 9 × 10 −4 K −1 at 134 K to 2 0 × 10 −4 K −1 at 350 K. In the transition region between 350 K and 405 K, it increases drastically before reaching the rubbery state, at an average value of 5 8 × 10 −4 K −1 . In Figure 6, linear fittings are then executed between 134 K and 350 K and between 405 K and 485 K. The experimental T g obtained at the intersection is thus 380 K.
The simulation/experiment match is good. The experimental and calculated CVTEs are similar in the glassy state (1 − 2 × 10 −4 and 2 × 10 −4 K −1 ) and in the good range in the rubbery state (5 8 × 10 −4 and 8 × 10 −4 K −1 ). The temperature domain for the transition state is much larger in the simulations (175 K) than in the experiments (55 K), likely because this polymer is not suitable for MD at temperatures above 500 K. This adds uncertainty to the determination of the simulated glass transition. Nevertheless, the simulated T g (390 ± 8 K) compares very well with both experimental glass transition temperatures of 391 ± 1 K measured by differential scanning calorimetry and through TMA (380 K). This is a good result considering the t − T dependence of the results. Actually, simulations are performed at cooling rates (q sim = 4 17 × 10 8 ) that are 10 8 orders of magnitude above that of TMA (q exp,TMA = 5 deg/min) and DSC (q exp,DSC = 10 deg/min). A rearrangement of the Williams-Landel-Ferry (WLF) equation (3), based on the t − T superposition principle, can express the expected difference between simulated and experimental T g [46]: where C1 and C2 are the WLF equation parameters. Using "universal" values (C1 = 17 44 K and C2 = 51 6 K), the simulated T g should be found 42.9 K above the TMA one and 40.0 K above the DSC one. This is above what we determine, and we may have to find more appropriate WLF parameters. These parameters depend on the fraction of free volume within the polymer and on the CVTE of this free volume, and they may be related to the polymer density. Since the matching between experimental and simulated densities is not perfect, there may be some level of porosity within the experimental samples that slightly modifies the CVTE vs. T behavior and then modifies the determination of T g,exp . Conversely, C1 and C2 values are probably not correct for our system which is highly cross-linked. For instance, in the DGEBA-TETA system, if the cross-linking density is high, then the free volume fraction is low (1.6-1.8%) and C1 equals 24-27 [47]. With a lower functionality and a longer aliphatic chain, the more flexible DGEBA-diaminomethylpentane system exhibits a higher fraction of free volume (4.3%) and a lower C1 value (10) [48]. ΔT g would then be much higher in the latter example. Thus, C1 and C2 should be better determined for the DGEBA-EDA system with a series of experiments and/or calculations at different cooling rates, as it has been done for other DGEBA-primary amine systems [49].

Conclusions
In our effort to model epoxy polymers, we have created a bulk polymer model starting from a stoichiometric liquid mixture of the two reactants DGEBA and EDA, EDA being a small liquid diamine molecule with a functionality of 4, similar enough to amines used in the industry, such as DETA and TETA. Six model structures were created through a stepwise procedure. We developed a cross-linking code that identifies any close reactive centers (reaction radius 3 Å) and creates bonds. The global polymerization procedure comprises also all-atom molecular dynamics (GAFF 1.8 force field) equilibration steps alternating with cross-linking. For the first five models, we performed 1 ns-long 700 K MD simulations between two cross-linking steps. In the latter steps, the diffusion of the monomers in the partially reticulated polymer could limit the cross-linking of the remaining reactive sites.
To favor the diffusion of the monomers, we increased the temperature of the MD simulations to 900 K, but there was no significant increase in the conversion degree that is finally higher than 90% for the six boxes. Thanks to additional molecular dynamics simulations, we were furthermore able to equilibrate and calculate the properties for the liquid mixture, as well as for those of the intermediate semipolymerized phases and the final polymer model. The properties (the density and glass transition temperature) of the final polymer model compare very well with experimental values, despite the uncertainty of the fitting of the rubbery state temperature region. Finally, the proposed bulk polyepoxy 3D-network is a good model, in agreement with its experimental counterpart showing the validity of the cross-linking procedure. Work is in progress to propose a surface model that will be useful to study the surface reactivity of epoxy polymers.

Data Availability
The data related to the present work are available on demand to the corresponding author.

Conflicts of Interest
The authors declare that they have no conflicts of interest.