polymers Article Characterization of Cure Behavior in Epoxy Using Molecular Dynamics Simulation Compared with Dielectric Analysis and DSC Shuang Yan * , Wolfgang Verestek , Harald Zeizinger and Siegfried Schmauder ���������� ������� Citation: Yan, S.; Verestek, W.; Zeizinger, H.; Schmauder, S. Characterization of Cure Behavior in Epoxy Using Molecular Dynamics Simulation Compared with Dielectric Analysis and DSC. Polymers 2021, 13, 3085. https://doi.org/10.3390/ polym13183085 Academic Editor: Keon-Soo Jang Received: 8 August 2021 Accepted: 10 September 2021 Published: 13 September 2021 Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affil- iations. Copyright: © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/). Institute for Materials Testing, Materials Science and Strength of Materials, University of Stuttgart, Pfaffenwaldring 32, 70569 Stuttgart, Germany; wolfgang.verestek@imwf.uni-stuttgart.de (W.V.); zeizinger@arcor.de (H.Z.); siegfried.schmauder@imwf.uni-stuttgart.de (S.S.) * Correspondence: shuang.yan@imwf.uni-stuttgart.de Abstract: The curing behavior of a thermosetting material that influences the properties of the mate- rial is a key issue for predicting the changes in material properties during processing. An empirical equation can describe the reaction kinetics of the curing behavior of an investigated material, which is usually estimated using experimental methods. In this study, the curing process of an epoxy resin, the polymer matrix in an epoxy molding compound, is computed concerning thermal influence using molecular dynamics. Furthermore, the accelerated reaction kinetics, which are influenced by an increased reaction cutoff distance, are investigated. As a result, the simulated crosslink density with various cutoff distances increases to plateau at a crosslink density of approx. 90% for the investigated temperatures during curing time. The reaction kinetics are derived according to the numerical results and compared with the results using experimental methods (dielectric analysis and differential scanning calorimetry), whereby the comparison shows a good agreement between experiment and simulation. Keywords: epoxy; curing; reaction kinetics; molecular dynamics; dielectric analysis (DEA); differen- tial scanning calorimetry (DSC) 1. Introduction Thermoset materials are widely used in the industrial sector because of their excellent mechanical properties at high temperatures and good chemical resistance. The non-cured thermosetting resin contains many monomers, which crosslink from a specific temperature and build three-dimensional macromolecules during processing. Since it causes changes in dynamical viscosity and temperature in the material, it is essential to know how the cure behavior changes during processing, e.g., the injection molding process. Several experi- mental methods are available for monitoring and investigating the thermoset material’s cure behavior: differential scanning calorimetry (DSC), dielectric analysis (DEA), Fourier transform infrared spectroscopy (FTIR), near-infrared spectroscopy (NIRs), and ultra-sonic techniques. However, using the measurement technology to monitor curing involves a high cost and effort and, in terms of the manufacturing process, the sensors cannot be adopted at every position where the curing state in the material needs to be measured. The molecular changes in morphological structure caused by the curing, which influence the material properties, are more easily predicted and analyzed with the help of molecular dynamics (MD) simulation. Molecular dynamics simulation is commonly used in thermoset polymers for ma- terial design and to predict material properties. The previous studies [1–15] focused on predicting the thermomechanical properties (e.g., glass transition temperature, Young’s modulus, thermal conductivity, and compressibility) in already cured molecular struc- tures of the epoxy polymers using MD simulation. The molecular structures (Li et al. [4]) Polymers 2021, 13, 3085. https://doi.org/10.3390/polym13183085 https://www.mdpi.com/journal/polymers https://www.mdpi.com/journal/polymers https://www.mdpi.com https://orcid.org/0000-0002-7869-8934 https://orcid.org/0000-0002-5155-0630 https://orcid.org/0000-0001-6022-3853 https://doi.org/10.3390/polym13183085 https://doi.org/10.3390/polym13183085 https://creativecommons.org/ https://creativecommons.org/licenses/by/4.0/ https://creativecommons.org/licenses/by/4.0/ https://doi.org/10.3390/polym13183085 https://www.mdpi.com/journal/polymers https://www.mdpi.com/article/10.3390/polym13183085?type=check_update&version=2 Polymers 2021, 13, 3085 2 of 16 and the crosslink density (Schichtel et al. [10]) strongly influence the thermomechanical properties. Okabe et al. [9] investigated the curing behavior of diglycidyl ether bisphenol A (DGEBA), tetra glycidyl diamino diphenyl methane (TGDDM), and eight polymers of diglycidyl ether bisphenol A (DGEBA8) as base resins, with 4′4-diaminodiphenylsulphone (44DDS) and 3′3-diaminodiphenylsulphone (33DDS) used as curing agents. The simu- lations showed an acceptable comparability of curing curves in different resin mixtures computed using MD simulation with experimental DSC measurements. Although the curing behavior was simulated in this study, it focused more on the resulting mechanical properties caused by different functional groups. In the study from Li et al. [16], the model was extended to a hybrid molecular dynamics/Monte Carlo (MD/MC) approach with reaction probabilities based on reaction path energies. Reactions of phenol resin with formaldehyde are investigated, and several thermal-mechanical properties agree with the experimental values. Unger et al. [17] attempted to validate the molecular modeling of a pure epoxy resin (base resins: bisphenol A diglycidyl ether, 1,4-Butanediol diglycidyl ether, Alkyl (C12-C14) glycidyl ether; hardener: α,ω-poly(oxypropylene)diamine, 3-(aminomethyl)- 3,5,5-trimethyl-cyclohexanamine, 2,4,6-tris[(dimethylamino)methyl]phenol) according to the spectroscopic changes in correlation with crosslink density by means of in situ NIR, which showed a good agreement between numerical and experimental results. The above-mentioned studies have shown that MD can simulate the curing of epoxy networks. However, the cure behavior (reaction kinetics) for different material systems was not investigated in detail because the curing was considered simply as a modeling step, and only the thermomechanical properties were checked. Either that or united atom force fields, e.g., the DREIDING force field, were used, which has been shown by Li et al. [18] to result in different morphologies than all-atom force fields. Until now, a validation of molecular modeling for curing using experimental methods is still challenging due to the difference in time scale between the molecular dynamics simulation in the nanoseconds range and the experimental study in the seconds range. This work focuses on analyzing the cure behavior of an epoxy resin (base resin: bisphenol A diglycidyl ether (BADGE) and hardener: ethylenediamine (EN)) with consideration of the thermal influence of using MD simulation. In addition, the accelerated reaction kinetics influenced by an increased reaction cutoff distance are investigated. Finally, the numerical results are compared with experimental results derived by DSC and DEA, whereby DSC is only applicable in a laboratory environment and DEA is suitable for in situ cure monitoring. A commercial epoxy composite applied in industrial application is measured in our experimental investigations, whose polymer matrix (base resin and hardener) is modeled in the MD simulation. Reaction models as part of reaction kinetics for describing the cure mechanism are estimated simultaneously using MD, DSC and DEA methods for evaluating their comparability. 2. Materials and Methods 2.1. Material The resin system, also referred to as EP in the following text, consists of epoxy resin bisphenol A diglycidyl ether (BADGE) and a hardener, ethylenediamine (EN), that crosslink via a polymerization reaction (see Table 1). In this molecular modeling, three possibilities for the crosslinking reaction baetween the functional groups are considered: (a) Reaction between the epoxide group (R1-CCO) and the primary amine group (R2-NH2); (b) Reaction between the epoxide group (R1-CCO) and the secondary amine group (R3- NH-R4); (c) Reaction between the epoxide group (R1-CCO) and the hydroxyl group (R5-OH). Polymers 2021, 13, 3085 3 of 16 Table 1. Overview of the resin system. Epoxy Resin Hardener Chemical formula bisphenol A diglycidyl ether (BADGE) ethylenediamine (EN) Molar mass 340.419 g/mol 60.1 g/mol Molecular formula C21H24O4 C2H8N2 Chemical structure Polymers 2021, 13, x FOR PEER REVIEW 3 of 17 The detailed sequence of the crosslinking reaction with molecular structures as ex- amples, is shown in Figure A1 in Appendix A. Table 1. Overview of the resin system. Epoxy Resin Hardener Chemical formula bisphenol A diglycidyl ether (BADGE) ethylenediamine (EN) Molar mass 340.419 g/mol 60.1 g/mol Molecular formula C21H24O4 C2H8N2 Chemical structure 2.2. Molecular Modeling For computing the chemical reactions (bond breaking, bond formation, and reactiv- ity) and predicting the material properties, there are two types of force fields: fixed force fields (also called classical force fields) and reactive force fields. For classical force fields (e.g., CVFF, PCFF, DREIDING, COMPASS), bonds are defined explicitly, while maintain- ing a low computing time and good accuracy [1,3,11,12,19–22]. In studies [8,23], a reactive force field (ReaxFF) is used to simulate the chemical reaction. The advantage of ReaxFF is that the bond formation is defined implicitly and happens dynamically during compu- ting, which needs a very high computing capacity and a high complexity of parameteri- zation. In this study, the crosslinking reactions between the functional groups are imple- mented in molecular modeling using a fixed force field (PCFF). The software EMC (Enhanced Monte Carlo), developed by Pieter J. in 't Veld [24], and LAMMPS Simulator (Large-scale Atomic/Molecular Massively Parallel, Sandia Na- tional Laboratories, USA) [25–27] were applied to prepare the initial molecular structures and molecular modeling. • Step 1. Generation of the initial molecular network of the uncured resin system. The first step was to generate the initial stochastic mixture of the uncured resin sys- tem using EMC software. EMC followed a combined Monte Carlo/MD approach. The fixed force field PCFF (Sun et al. [28]) was used to parameterize intra- and inter-molecular interactions between the atoms. The ratio of epoxy resin BADGE to the hardener EN was 2:1. The epoxy resin BADGE amount is NBADGE = 1817, which included the amount of the epoxide group N-CCO = 3634. The amount of the hardener EN NEN = 908, and includes the primary amine group N-NH2 = 1818. The initial density after creation by EMC was 1 g/cm³. As an initial equilibration, a short simulation sequence in microcanonical ensemble NVE (with limited maximum displacement per step), canonical ensemble NVT, and iso- thermal-isobaric ensemble NPT was performed in LAMMPS with a time step of 0.25fs, at a target temperature of 150 °C, and for a total time of 25ps. • Step 2. Thermodynamical equilibration of the initial uncured molecular network. The initial uncured molecular mixture was equilibrated for 150ps under periodic boundary conditions using NPT to achieve thermodynamical equilibrium. A time step size of 1.5fs was applied for these and all following simulations by utilizing the rRESPA algorithm (Plimpton et al. [26]) implemented in LAMMPS. The conditions (pressure P: 423 atm; temperature T: 150 °C, 170 °C, and 190 °C), which were applied as process parameters in DEA measurements in an actual injection mold, were implemented for the equilibration procedure. Table 2 lists the density of the uncured molecular system in a simulation box (including its dimension and volume) after equilibration at three temperatures. Polymers 2021, 13, x FOR PEER REVIEW 3 of 17 The detailed sequence of the crosslinking reaction with molecular structures as ex- amples, is shown in Figure A1 in Appendix A. Table 1. Overview of the resin system. Epoxy Resin Hardener Chemical formula bisphenol A diglycidyl ether (BADGE) ethylenediamine (EN) Molar mass 340.419 g/mol 60.1 g/mol Molecular formula C21H24O4 C2H8N2 Chemical structure 2.2. Molecular Modeling For computing the chemical reactions (bond breaking, bond formation, and reactiv- ity) and predicting the material properties, there are two types of force fields: fixed force fields (also called classical force fields) and reactive force fields. For classical force fields (e.g., CVFF, PCFF, DREIDING, COMPASS), bonds are defined explicitly, while maintain- ing a low computing time and good accuracy [1,3,11,12,19–22]. In studies [8,23], a reactive force field (ReaxFF) is used to simulate the chemical reaction. The advantage of ReaxFF is that the bond formation is defined implicitly and happens dynamically during compu- ting, which needs a very high computing capacity and a high complexity of parameteri- zation. In this study, the crosslinking reactions between the functional groups are imple- mented in molecular modeling using a fixed force field (PCFF). The software EMC (Enhanced Monte Carlo), developed by Pieter J. in 't Veld [24], and LAMMPS Simulator (Large-scale Atomic/Molecular Massively Parallel, Sandia Na- tional Laboratories, USA) [25–27] were applied to prepare the initial molecular structures and molecular modeling. • Step 1. Generation of the initial molecular network of the uncured resin system. The first step was to generate the initial stochastic mixture of the uncured resin sys- tem using EMC software. EMC followed a combined Monte Carlo/MD approach. The fixed force field PCFF (Sun et al. [28]) was used to parameterize intra- and inter-molecular interactions between the atoms. The ratio of epoxy resin BADGE to the hardener EN was 2:1. The epoxy resin BADGE amount is NBADGE = 1817, which included the amount of the epoxide group N-CCO = 3634. The amount of the hardener EN NEN = 908, and includes the primary amine group N-NH2 = 1818. The initial density after creation by EMC was 1 g/cm³. As an initial equilibration, a short simulation sequence in microcanonical ensemble NVE (with limited maximum displacement per step), canonical ensemble NVT, and iso- thermal-isobaric ensemble NPT was performed in LAMMPS with a time step of 0.25fs, at a target temperature of 150 °C, and for a total time of 25ps. • Step 2. Thermodynamical equilibration of the initial uncured molecular network. The initial uncured molecular mixture was equilibrated for 150ps under periodic boundary conditions using NPT to achieve thermodynamical equilibrium. A time step size of 1.5fs was applied for these and all following simulations by utilizing the rRESPA algorithm (Plimpton et al. [26]) implemented in LAMMPS. The conditions (pressure P: 423 atm; temperature T: 150 °C, 170 °C, and 190 °C), which were applied as process parameters in DEA measurements in an actual injection mold, were implemented for the equilibration procedure. Table 2 lists the density of the uncured molecular system in a simulation box (including its dimension and volume) after equilibration at three temperatures. The detailed sequence of the crosslinking reaction with molecular structures as exam- ples, is shown in Figure A1 in Appendix A. 2.2. Molecular Modeling For computing the chemical reactions (bond breaking, bond formation, and reactivity) and predicting the material properties, there are two types of force fields: fixed force fields (also called classical force fields) and reactive force fields. For classical force fields (e.g., CVFF, PCFF, DREIDING, COMPASS), bonds are defined explicitly, while maintaining a low computing time and good accuracy [1,3,11,12,19–22]. In studies [8,23], a reactive force field (ReaxFF) is used to simulate the chemical reaction. The advantage of ReaxFF is that the bond formation is defined implicitly and happens dynamically during computing, which needs a very high computing capacity and a high complexity of parameterization. In this study, the crosslinking reactions between the functional groups are implemented in molecular modeling using a fixed force field (PCFF). The software EMC (Enhanced Monte Carlo), developed by Pieter J. in ’t Veld [24], and LAMMPS Simulator (Large-scale Atomic/Molecular Massively Parallel, Sandia National Laboratories, USA) [25–27] were applied to prepare the initial molecular structures and molecular modeling. • Step 1. Generation of the initial molecular network of the uncured resin system. The first step was to generate the initial stochastic mixture of the uncured resin system using EMC software. EMC followed a combined Monte Carlo/MD approach. The fixed force field PCFF (Sun et al. [28]) was used to parameterize intra- and inter-molecular interactions between the atoms. The ratio of epoxy resin BADGE to the hardener EN was 2:1. The epoxy resin BADGE amount is NBADGE = 1817, which included the amount of the epoxide group N-CCO = 3634. The amount of the hardener EN NEN = 908, and includes the primary amine group N-NH2 = 1818. The initial density after creation by EMC was 1 g/cm3. As an initial equilibration, a short simulation sequence in microcanonical ensem- ble NVE (with limited maximum displacement per step), canonical ensemble NVT, and isothermal-isobaric ensemble NPT was performed in LAMMPS with a time step of 0.25 fs, at a target temperature of 150 ◦C, and for a total time of 25 ps. • Step 2. Thermodynamical equilibration of the initial uncured molecular network. The initial uncured molecular mixture was equilibrated for 150ps under periodic boundary conditions using NPT to achieve thermodynamical equilibrium. A time step size of 1.5 fs was applied for these and all following simulations by utilizing the rRESPA algo- rithm (Plimpton et al. [26]) implemented in LAMMPS. The conditions (pressure P: 423 atm; temperature T: 150 ◦C, 170 ◦C, and 190 ◦C), which were applied as process parameters in DEA measurements in an actual injection mold, were implemented for the equilibration procedure. Table 2 lists the density of the uncured molecular system in a simulation box (including its dimension and volume) after equilibration at three temperatures. Polymers 2021, 13, 3085 4 of 16 Table 2. The density of the initial molecular mixture in a simulation cell equilibrated at 150 ◦C, 170 ◦C, and 190 ◦C. T [◦C] Initial Density [g/cm3] Initial Dimension [nm3] Initial Volume [nm3] 150 0.988 10.4 × 10.4 × 10.4 1131.3 170 0.977 10.4 × 10.4 × 10.4 1144.2 190 0.962 10.5 × 10.5 × 10.5 1162.1 • Step 3. Identification of the reactive bonding and definition of criteria for the beginning of crosslinking. For crosslinking resin and hardener, the REACTER protocol [29,30] as implemented in LAMMPS was used. The cutoff distance Dcut, which correlated to the bond energy, was defined as a criterion for bond breaking and creation. Bond dissociation energy described the amount of energy that was required to split the bond homolytically. Numerically, the reactivity of the crosslinking reaction was changed through varying the cutoff distance Dcut. By increasing cutoff distance, the reactivity increased, and the duration until the material was fully crosslinked reduced. For example, during curing between epoxide and amine groups, if the distance between functional groups satisfied the boundary condition (≤cutoff distance Dcut), the bond C-O in the epoxide group will break up and link up to the bond N-H, which also breaks simultaneously. As a reference, the bond dissociation length of the newly created bonds (O-H, C-N, and C-O) was calculated using the class2 force field (see Table 3). The bond potentials Vbond between the atom i and j were generated with Equation (1) using class2 bond style presented by Sun [31]: Vbond,ij ( rij ) = K2· ( rij − rij.0 )2 + K3· ( rij − rij.0 )3 + K4· ( rij − rij.0 )4 (1) Table 3. The bond dissociation length of the newly created bonds. Bond-Type (i– j) Unit O–H C–N C–O Bond dissociation energy Vbond,ij [32] [kJ/mol] 463 305 358 Bond dissociation energy Vbond,ij [32] [kcal/mol] 110.7 72.9 85.56 rij.0 [Å] 0.965 1.452 1.42 K2 - 532.5062 327.1657 400.3954 K3 - −1282.905 −547.899 −835.1951 K4 - 2004.7658 526.5 1313.0142 Bond dissociation length rij [Å] 1.4858 2.0717 1.9373 Parameter rij.0 presents the equilibrium bond distance. K2, K3, and K4 are bond coefficients. This molecular modeling of crosslinking performed using the classic force field showed several advantages. It needed less computing capacity and performed more efficiently. The temperature influence and the effect of cutoff distance Dcut on the cure ki- netics were analyzed. The curing process was simulated at temperatures of 150 ◦C, 170 ◦C, and 190 ◦C while varying the cutoff distance Dcut of the O-H bond from 1.5 Å→2.5 Å. 2.3. Dielectric Analysis The dielectric analysis is a measuring procedure to study the material’s response to an applied electric field. The electrical current Ielectr(f) flowing through an electric capacitor, in which the epoxy material as dielectric material was located, was measured as a response to an alternating electric field depending on field frequency f. This analyzed the changes in dielectric properties, e.g., the relative permittivity εr of the polymer, depending on the changes in molecular and morphological structures in the material. Therefore, the the dialectric behavior of thermosets can be used to monitor the curing process which causes changes in dielectric properties. Permittivity ε, determined by the material polarization in Polymers 2021, 13, 3085 5 of 16 response to the electrical field, is a physical quantity describing the interaction between the electric field and the dielectric material. Relative permittivity εr is defined as the ratio of dielectric permittivity ε to vacuum permittivity ε0. The electric field polarized the dipoles and freed the charges in the investigated material, leading to a change in relative permittivity related to the curing. The relative permittivity was calculated according to the measured electrical current and electrical voltage. The electrical current Ielectr resulted from dipole interaction in the capacitor’s electric field and the mobility of ionized molecules. The relative permittivity εr of the investigated material was a complex quantity with angular frequency ω = 2πf, shown in Equation (2). εr(ω) = ε′(ω)− i·ε′′ (ω) (2) The term ε’ is the real part of relative permittivity that represents stored electric energy in the material. The imaginary part ε” of the relative permittivity (dielectric loss) describes the lost energy of the applied external electric field. In this study, dielectric measurements were performed using a sensor TMC 16/3 (NETZSCH-Gerätebau GmbH, Selb, Germany) with 1 kHz at a pressure of 423 atm and temperatures of 150 ◦C, 170 ◦C, and 190 ◦C. Figure 1 shows an example of the evaluation procedure of a DEA measuring result obtained at the temperature of 150 ◦C. The change of relative permittivity εr during heating is shown in Figure 1a. Firstly the relative permittivity εr increases until t0, which indicates the melting of the materials and the increased amount of available charge carriers in the material. After the onset of curing, the polymer system transits from a highly viscous to solid-state, and the content of free movable charges reduces, leading to a decrease in the relative permittivity εr and an increase in crosslink density, meaning that εr ~ 1/α. The maximal relative permittivity εr,0 at t0 means the curing begins (α ≈ 0%). After the material is almost fully cured, the relative permittivity εr reaches a plateau. Termination of curing (α ≈ 100%) is determined at time t1 when the relative permittivity’s change rate is less than 1%. Thereby, the value of relative permittivity εr,1 is determined at time t1. Based on the relative permittivity’s change in correlation with crosslink density α during curing, the crosslink density is estimated as a function of heating time according to Equation (3), shown in Figure 1b. In the further step, the reaction rate dα/dt in dependency on the crosslink density α, shown as dots in Figure 1c, can be derived regarding the results of the last step to determine the unknown parameter in reaction model f(α). Whereby the black curve of f(α) is the fitted line regarding the values of dα/dt vs. α. αDEA = (εr − εr,∆)/εr,∆ εr,∆ = εr,1 − εr,0 (3) αDEA is the crosslink density determined using DEA. The term εr,0 means the maximal relative permittivity occurring at t0; εr,1 is the relative permittivity determined at time t1, when the reaction is terminated, and εr,∆ is the value of the difference between εr,0 and εr,1. 2.4. Differential Scanning Calorimetry DSC is a classic method for analyzing the correlation between the curing amount, temperature, and time, based on the amount of heat released from crosslinking reactions. The equipment DSC Q100 (TA Instruments–Waters LLC, New Castle, DE, USA) was used for studying the reaction kinetics. Non-isothermal measurements were performed with heating rates of 2 ◦C/min, 5 ◦C/min, 10 ◦C/min, 15 ◦C/min, and 20 ◦C/min with nitrogen N2 as the purge gas. Polymers 2021, 13, 3085 6 of 16 Polymers 2021, 13, x FOR PEER REVIEW 6 of 17 Figure 1. Evaluation procedure of DEA result: (a) maximum and minimum of relative permittivity, (b) crosslink density based on DEA result, and (c) reaction rate in dependency of crosslink density. 2.4. Differential Scanning Calorimetry DSC is a classic method for analyzing the correlation between the curing amount, temperature, and time, based on the amount of heat released from crosslinking reactions. The equipment DSC Q100 (TA Instruments–Waters LLC, New Castle, DE, USA) was used for studying the reaction kinetics. Non-isothermal measurements were performed with heating rates of 2 °C/min, 5 °C/min, 10 °C/min, 15 °C/min, and 20 °C/min with nitrogen N2 as the purge gas. The thermally stimulated depolarization currents method [33,34] and DSC are appli- cable for an experimental determination of the activation energy, while the DSC method has been employed in our study. Concerning the results of DSC measurements, the acti- vation energy EA and the pre-exponential A could be determined using the Kissinger method [35]. For predicting the unknown parameter in reaction model f(α), the following procedure is performed. During the curing process, the change of heat flow in depend- ency of time and temperature is detected by DSC shown in Figure 2a. According to the resultant change of heat amount, temperature- and time-dependent crosslink density is a ratio of the released heat amount Hrelease to the maximum heat amount Hmax in a fully cross- linked state, using Equation (4) (see Figure 2b). Figure 2c shows the reaction rate in de- pendency on the crosslinking degree to estimate the unknown parameter in the reaction model f(α): 𝛼 = 𝐻 𝐻⁄ (4) The term αDSC represents the crosslink density determined using DSC. The Hrelease is the released heat amount by the material in the partially cured state, and Hmax is the max- imum heat amount released by the material in the fully cured state. Figure 1. Evaluation procedure of DEA result: (a) maximum and minimum of relative permittivity, (b) crosslink density based on DEA result, and (c) reaction rate in dependency of crosslink density. The thermally stimulated depolarization currents method [33,34] and DSC are applica- ble for an experimental determination of the activation energy, while the DSC method has been employed in our study. Concerning the results of DSC measurements, the activation energy EA and the pre-exponential A could be determined using the Kissinger method [35]. For predicting the unknown parameter in reaction model f(α), the following procedure is performed. During the curing process, the change of heat flow in dependency of time and temperature is detected by DSC shown in Figure 2a. According to the resultant change of heat amount, temperature- and time-dependent crosslink density is a ratio of the released heat amount Hrelease to the maximum heat amount Hmax in a fully crosslinked state, using Equation (4) (see Figure 2b). Figure 2c shows the reaction rate in dependency on the crosslinking degree to estimate the unknown parameter in the reaction model f(α): αDSC = Hrelease/Hmax (4) The term αDSC represents the crosslink density determined using DSC. The Hrelease is the released heat amount by the material in the partially cured state, and Hmax is the maximum heat amount released by the material in the fully cured state. Polymers 2021, 13, 3085 7 of 16 Polymers 2021, 13, x FOR PEER REVIEW 7 of 17 Figure 2. Evaluation procedure of DSC result measured under the non-isothermal condition: (a) change of heat flow dur- ing the curing process with heating rate 10 °C/min, (b) crosslink density, and (c) reaction rate in dependency of curing amount. 3. Results and Discussion 3.1. Curing Behavior Once the initial uncured molecular mixtures were finished to equilibrate with NPT, the crosslinking simulations were performed with varied cutoff distances for O-H bonds from 1.48 Å up to 2.5 Å at 150 °C, 170 °C, and 190 °C (shown in Figure 3a–c). The number of completed reactions is counted during simulation. Since the epoxide group -CCO par- ticipates in each of the modeled crosslinking reactions, Equation (5) is derived for estimat- ing the crosslink density as the ratio of the number of already utilized epoxide groups Nu,CCO to the total number of epoxide groups Ntotal,CCO. The computing time and effort, with a calculated bond dissociation length (Dcut = 1.48 Å), are enormous until the material is fully cured as obtained by the MD simulation. Then, further simulations are carried out with the increasing cutoff distance Dcut of 1.48 Å up to 2.5 Å, which increases the reactivity of the chemical reactions, respectively. Based on computing results, the cutoff distance has a significant influence on the cur- ing reaction’s activity. At three temperatures, the molecules crosslink very slowly with a low cutoff distance of 1.5 Å, and it is to be expected that a significant change in crosslink density is only visible over a very long simulation time. However, since the cutoff distance Figure 2. Evaluation procedure of DSC result measured under the non-isothermal condition: (a) change of heat flow during the curing process with heating rate 10 ◦C/min, (b) crosslink density, and (c) reaction rate in dependency of curing amount. 3. Results and Discussion 3.1. Curing Behavior Once the initial uncured molecular mixtures were finished to equilibrate with NPT, the crosslinking simulations were performed with varied cutoff distances for O-H bonds from 1.48 Å up to 2.5 Å at 150 ◦C, 170 ◦C, and 190 ◦C (shown in Figure 3a–c). The number of completed reactions is counted during simulation. Since the epoxide group -CCO participates in each of the modeled crosslinking reactions, Equation (5) is derived for estimating the crosslink density as the ratio of the number of already utilized epoxide groups Nu,CCO to the total number of epoxide groups Ntotal,CCO. The computing time and effort, with a calculated bond dissociation length (Dcut = 1.48 Å), are enormous until the material is fully cured as obtained by the MD simulation. Then, further simulations are carried out with the increasing cutoff distance Dcut of 1.48 Å up to 2.5 Å, which increases the reactivity of the chemical reactions, respectively. Based on computing results, the cutoff distance has a significant influence on the curing reaction’s activity. At three temperatures, the molecules crosslink very slowly with a low cutoff distance of 1.5 Å, and it is to be expected that a significant change in crosslink density is only visible over a very long simulation time. However, since the cutoff distance is set to larger than 1.6 Å, the crosslinking process becomes significantly faster. Additionally, Polymers 2021, 13, 3085 8 of 16 the maximal achieved crosslinking amount also increases. A similar effect of cutoff distance on the crosslinking behavior is also observed in the report by Unger et al. [17]. Polymers 2021, 13, x FOR PEER REVIEW 8 of 17 is set to larger than 1.6 Å, the crosslinking process becomes significantly faster. Addition- ally, the maximal achieved crosslinking amount also increases. A similar effect of cutoff distance on the crosslinking behavior is also observed in the report by Unger et al. [17]. (a) (b) (c) Figure 3. Curves of crosslink density depending on curing time concerning the influences of cutoff distances (1.48 Å→2.5 Å) at (a) 150 °C, (b) 170 °C, and (c) 190 °C. Generally, a higher environment temperature (e.g., 190 °C) excites stronger molecu- lar oscillations, which cause a higher molecular diffusion in the simulation system. Thus, while the molecular movement becomes more active at a higher temperature of 190 °C than at a lower temperature of 150 °C, the epoxy resin, over time, has a higher probability for meeting a hardener and linking up. Respectively, an achievement of a higher crosslink density at a higher temperature is expected, which is demonstrated in our simulation study with the cutoff distance of 1.6 Å. However, the influence of the rising environmental temperature on the crosslink density is no more significantly recognized in respect of nu- merical results computed with a cutoff distance of ≥ 1.7 Å. As described in step 2 in Section 2.2, the molecular mixture is previously equilibrated at the given temperature, and all molecules distribute homogeneously in space. Since the diffusion of molecules is a long-term process which only occurs at a lower cutoff distance (e.g., 1.6 Å), indicating a longer curing time until the maximum crosslink density is achieved, the influence of the increasing temperature (150 °C→190 °C) on diffusion and the corresponding crosslink density is observed in Figure 3. By increasing the cutoff dis- tance (≥ 1.7 Å), the molecules have a higher probability of reacting with the neighboring reaction partners in a homogenized polymer mixture, while the simulated curing time is reduced, therefore preventing molecular diffusion. As a result, excited intermolecular in- teractions are not significantly observable due to the increasing temperature, which leads to a similar crosslink density value simulated with the equal cutoff distance at varying temperatures (150 °C→190 °C). Additionally, the molecular simulation demonstrates that a fully cured molecular network (α ≈ 100%) is hardly realizable. During the curing process, more crosslinked mo- lecular networks are created and grow continuously to a larger size. The large molecules restrict the remaining epoxy resin and hardener from free movement and finding the re- action partner. 𝛼 = 𝑁 ,𝑁 . (5) The term αMD represents the crosslink density determined using molecular dynamics simulation. Nu,CCO is the number of epoxide groups already utilized during curing, and Ntotal,CCO is the total number of epoxide groups. Table 4 lists the densities of the molecular mixtures at the end of the curing simula- tions for cutoff distances (1.48 Å→2.5 Å) at 150 °C, 170 °C, and 190 °C. Referenced to the Figure 3. Curves of crosslink density depending on curing time concerning the influences of cutoff distances (1.48 Å→2.5 Å) at (a) 150 ◦C, (b) 170 ◦C, and (c) 190 ◦C. Generally, a higher environment temperature (e.g., 190 ◦C) excites stronger molecular oscillations, which cause a higher molecular diffusion in the simulation system. Thus, while the molecular movement becomes more active at a higher temperature of 190 ◦C than at a lower temperature of 150 ◦C, the epoxy resin, over time, has a higher probability for meeting a hardener and linking up. Respectively, an achievement of a higher crosslink density at a higher temperature is expected, which is demonstrated in our simulation study with the cutoff distance of 1.6 Å. However, the influence of the rising environmental temperature on the crosslink density is no more significantly recognized in respect of numerical results computed with a cutoff distance of ≥ 1.7 Å. As described in step 2 in Section 2.2, the molecular mixture is previously equilibrated at the given temperature, and all molecules distribute homogeneously in space. Since the diffusion of molecules is a long-term process which only occurs at a lower cutoff distance (e.g., 1.6 Å), indicating a longer curing time until the maximum crosslink density is achieved, the influence of the increasing temperature (150 ◦C→190 ◦C) on diffusion and the corresponding crosslink density is observed in Figure 3. By increasing the cutoff distance (≥ 1.7 Å), the molecules have a higher probability of reacting with the neighboring reaction partners in a homogenized polymer mixture, while the simulated curing time is reduced, therefore preventing molecular diffusion. As a result, excited intermolecular interactions are not significantly observable due to the increasing temperature, which leads to a similar crosslink density value simulated with the equal cutoff distance at varying temperatures (150 ◦C→190 ◦C). Additionally, the molecular simulation demonstrates that a fully cured molecular network (α ≈ 100%) is hardly realizable. During the curing process, more crosslinked molecular networks are created and grow continuously to a larger size. The large molecules restrict the remaining epoxy resin and hardener from free movement and finding the reaction partner. αMD = Nu,CCO Ntotal.CCO (5) The term αMD represents the crosslink density determined using molecular dynamics simulation. Nu,CCO is the number of epoxide groups already utilized during curing, and Ntotal,CCO is the total number of epoxide groups. Table 4 lists the densities of the molecular mixtures at the end of the curing simulations for cutoff distances (1.48 Å→2.5 Å) at 150 ◦C, 170 ◦C, and 190 ◦C. Referenced to the initial densities of the uncured molecular system in Table 2, the molecular mixture’s density at the end of simulation increases when the achieved crosslink density grows by increasing the cutoff distances shown in Figure 3. Additionally, the density decreases during the Polymers 2021, 13, 3085 9 of 16 increasing temperature of 150 ◦C→190 ◦C, which correlates to the thermal expansion of the molecular system influenced by a higher temperature. Table 4. Density of the molecular mixture at the end of the crosslinking simulation for Dcut = 1.48 Å→2.5 Å at 150 ◦C, 170 ◦C, and 190 ◦C. T [◦C] Density [g/cm3] Dcut = 1.48 Å Dcu = 1.5 Å Dcu = 1.6 Å Dcu = 1.7 Å Dcu = 1.8 Å Dcu = 2.5 Å 150 0.984 0.992 1.050 1.077 1.078 1.080 170 0.980 0.982 1.058 1.074 1.075 1.077 190 0.960 0.963 1.028 1.063 1.064 1.068 Changes in the molecular structures regarding the result with a cutoff distance of 1.8 Å at 190 ◦C are schematically shown as an example in Figure 4. The colors of electrically charged atoms change according to changes in the atoms’ charges during bond breaking and formation. Polymers 2021, 13, x FOR PEER REVIEW 9 of 17 initial densities of the uncured molecular system in Table 2, the molecular mixture’s den- sity at the end of simulation increases when the achieved crosslink density grows by in- creasing the cutoff distances shown in Figure 3. Additionally, the density decreases dur- ing the increasing temperature of 150 °C→190 °C, which correlates to the thermal expan- sion of the molecular system influenced by a higher temperature. Table 4. Density of the molecular mixture at the end of the crosslinking simulation for Dcut = 1.48 Å→2.5 Å at 150 °C, 170 °C, and 190 °C. T [°C] Density [g/cm³] Dcut = 1.48 Å Dcu = 1.5 Å Dcu = 1.6 Å Dcu = 1.7 Å Dcu = 1.8 Å Dcu = 2.5 Å 150 0.984 0.992 1.050 1.077 1.078 1.080 170 0.980 0.982 1.058 1.074 1.075 1.077 190 0.960 0.963 1.028 1.063 1.064 1.068 Changes in the molecular structures regarding the result with a cutoff distance of 1.8 Å at 190 °C are schematically shown as an example in Figure 4. The colors of electrically charged atoms change according to changes in the atoms’ charges during bond breaking and formation. (a) (b) Figure 4. Changes in molecular structures during the curing process with the cutoff distance of 1.8 Å at 190 °C: (a) total view of the simulation box with periodic boundary, and (b) growth of the molecule as an example. The reaction speed simulated by molecular dynamics is significantly higher than the experimentally determined value. In contrast to the experiment in which a larger test spec- imen is considered in millimeter scale, the material volume in nanometer-scale is taken into account in MD simulation, shown in Figure 5. An ideal thermodynamic state is set in molecular dynamics, under a constant pressure and temperature with a periodic bound- ary condition, meaning that there is no heat transfer within the material and no heat loss occurs between the edge and the environment. In the experimental analysis, heat transfer within the material and thermal exchange within the environment causes the imbalance of the temperature distribution across the layer thickness in the polymer matrix from EP. Figure 4. Changes in molecular structures during the curing process with the cutoff distance of 1.8 Å at 190 ◦C: (a) total view of the simulation box with periodic boundary, and (b) growth of the molecule as an example. The reaction speed simulated by molecular dynamics is significantly higher than the experimentally determined value. In contrast to the experiment in which a larger test specimen is considered in millimeter scale, the material volume in nanometer-scale is taken into account in MD simulation, shown in Figure 5. An ideal thermodynamic state is set in molecular dynamics, under a constant pressure and temperature with a periodic boundary condition, meaning that there is no heat transfer within the material and no heat loss occurs between the edge and the environment. In the experimental analysis, heat transfer within the material and thermal exchange within the environment causes the imbalance of the temperature distribution across the layer thickness in the polymer matrix from EP. Additionally, monomers are homogeneously distributed in the simulation volume, leading to a shorter balancing time of the molecular movements and a faster reaction. Furthermore, besides pure epoxy resin and hardener, the commercial EP material investigated in the Polymers 2021, 13, 3085 10 of 16 experiment contains additives such as the reaction inhibitor, which delays and slows down the curing reaction. Polymers 2021, 13, x FOR PEER REVIEW 10 of 17 Additionally, monomers are homogeneously distributed in the simulation volume, lead- ing to a shorter balancing time of the molecular movements and a faster reaction. Further- more, besides pure epoxy resin and hardener, the commercial EP material investigated in the experiment contains additives such as the reaction inhibitor, which delays and slows down the curing reaction. (a) (b) Figure 5. Macro- and micro-structure of the epoxy resin: (a) Nano-CT of the epoxy resin, and (b) view of simulation cell in molecular dynamics simulation. 3.2. Reaction Kinetics Using Molecular Dynamics Simulation Reaction kinetics describe the rate of the chemical reaction, which considers temper- ature, pressure and concentration dependency: k(T), f(α), and h(p) (see Equation (6)). 𝑑𝛼/𝑑𝑡 = 𝑘(𝑇) ∙ 𝑓(𝛼) ∙ ℎ(𝑝) (6) Supposing that the pressure p is kept constant during the chemical reaction process, the pressure dependency can be ignored, and Equation (6) is simplified into Equations (7) and (8). 𝑑𝛼/𝑑𝑡 = 𝑘(𝑇) ∙ 𝑓(𝛼) (7) The temperature dependency, k(T) = A·exp(-EA/(Rgas·T)), generally known as the Ar- rhenius equation, consists of the material-related activation energy EA and the pre-expo- nential term A: 𝑑𝛼/𝑑𝑡 = 𝐴 ∙ 𝑒𝑥𝑝(−𝐸 /(𝑅 ∙ 𝑇)) ∙ 𝑓(𝛼). (8) The reaction model f(α) describes the reaction mechanisms in dependency on the pro- portion of crosslinked molecules. There are already empirical equations available (such as nth-order catalytic, Avrami–Erofeev, Prout–Tompkins, and Sestak–Berggren reaction model [36,37]), whose unknown factors can be determined by regression methods. Since the empirical model (nth-order catalytic) best fits the numerical results for this investi- gated material, only this model is considered further in our study. Concerning the difference in the time scale between the molecular dynamics simula- tion in the nanoseconds range and the experimental study in the seconds range, the eval- Figure 5. Macro- and micro-structure of the epoxy resin: (a) Nano-CT of the epoxy resin, and (b) view of simulation cell in molecular dynamics simulation. 3.2. Reaction Kinetics Using Molecular Dynamics Simulation Reaction kinetics describe the rate of the chemical reaction, which considers tempera- ture, pressure and concentration dependency: k(T), f(α), and h(p) (see Equation (6)). dα/dt = k(T)· f (α)·h(p) (6) Supposing that the pressure p is kept constant during the chemical reaction process, the pressure dependency can be ignored, and Equation (6) is simplified into Equations (7) and (8). dα/dt = k(T)· f (α) (7) The temperature dependency, k(T) = A·exp(-EA/(Rgas·T)), generally known as the Arrhenius equation, consists of the material-related activation energy EA and the pre- exponential term A: dα/dt = A·exp ( −EA/ ( Rgas·T )) · f (α). (8) The reaction model f(α) describes the reaction mechanisms in dependency on the proportion of crosslinked molecules. There are already empirical equations available (such as nth-order catalytic, Avrami–Erofeev, Prout–Tompkins, and Sestak–Berggren reaction model [36,37]), whose unknown factors can be determined by regression methods. Since the empirical model (nth-order catalytic) best fits the numerical results for this investigated material, only this model is considered further in our study. Concerning the difference in the time scale between the molecular dynamics simu- lation in the nanoseconds range and the experimental study in the seconds range, the evaluation procedure for determining the reaction model f(α) regarding the experimental results can also be applied for the calculation based on the simulation results. As shown Polymers 2021, 13, 3085 11 of 16 in Figure 6a–c, the time dependent reaction rate is calculated based on crosslink density depending on curing time, which is simulated using various cutoff distances and tem- peratures in Figure 3. With a large cutoff (2.5 Å), the maximal reaction rate is achieved immediately after the beginning of curing. As the cutoff distance decreases, the maximal reaction rate decreases until no reaction occurs during the simulation time. The data of reaction rate and the associated crosslink density at the same curing time are transferred to Figure 6d–f to further explain the determination of the unknown parameters in reaction model f(α). After that, the reaction rate curve is normalized according to the maximum reaction rate, related to the associated cutoff distance and temperature. Since the tempera- ture dependency is considered by factor k(T) in Equation (7), the data (normalized dα/dt vs. α), simulated with the same cutoff distance at different temperatures, are merged for the subsequent regression procedure (see Figure 7). Polymers 2021, 13, x FOR PEER REVIEW 11 of 17 uation procedure for determining the reaction model f(α) regarding the experimental re- sults can also be applied for the calculation based on the simulation results. As shown in Figure 6a–c, the time dependent reaction rate is calculated based on crosslink density de- pending on curing time, which is simulated using various cutoff distances and tempera- tures in Figure 3. With a large cutoff (2.5 Å), the maximal reaction rate is achieved imme- diately after the beginning of curing. As the cutoff distance decreases, the maximal reac- tion rate decreases until no reaction occurs during the simulation time. The data of reac- tion rate and the associated crosslink density at the same curing time are transferred to Figure 6d–f to further explain the determination of the unknown parameters in reaction model f(α). After that, the reaction rate curve is normalized according to the maximum reaction rate, related to the associated cutoff distance and temperature. Since the temper- ature dependency is considered by factor k(T) in Equation (7), the data (normalized dα/dt vs. α), simulated with the same cutoff distance at different temperatures, are merged for the subsequent regression procedure (see Figure 7). (a) (b) (c) (d) (e) (f) Figure 6. (a–c) Curves of reaction rate depending on curing time and (d–f) curves of reaction rate depending on crosslink density concerning the influence of cutoff distance (1.48 Å→2.5 Å) at 150 °C, 170 °C, and 190 °C. Figure 7 shows the normalized reaction rate in the dependency of crosslink density computed using different cutoff distances. The plot of reaction rate scatters less by increas- ing cutoff distance. However, supposing that the cutoff distance is set too high, such as at 1.9 Å, the chemical reaction takes place too fast so that the maximum reaction rate is al- ready reached, even at a low crosslink density (Figure 7d). The reaction model f(α) is esti- mated based on the numerical results (normalized dα/dt vs. α) with the cutoff distance of 1.8 Å, whereby the regression result of the reaction model (nth-order catalytic) is pre- sented in Figure 7c. Figure 6. (a–c) Curves of reaction rate depending on curing time and (d–f) curves of reaction rate depending on crosslink density concerning the influence of cutoff distance (1.48 Å→2.5 Å) at 150 ◦C, 170 ◦C, and 190 ◦C. Figure 7 shows the normalized reaction rate in the dependency of crosslink density computed using different cutoff distances. The plot of reaction rate scatters less by increas- ing cutoff distance. However, supposing that the cutoff distance is set too high, such as at 1.9 Å, the chemical reaction takes place too fast so that the maximum reaction rate is already reached, even at a low crosslink density (Figure 7d). The reaction model f(α) is estimated based on the numerical results (normalized dα/dt vs. α) with the cutoff distance of 1.8 Å, whereby the regression result of the reaction model (nth-order catalytic) is presented in Figure 7c. 3.3. Comparison between Simulation and Experimental Results 3.3.1. Reaction Kinetics Using DEA DEA measurements were performed under the pressure of 423 atm and at temperatures of 150 ◦C, 170 ◦C, and 190 ◦C. The change in reaction rate depending on the crosslinking amount, which is calculated according to the relative permittivity’s changes at different temperatures, is plotted in Figure 8. At the beginning of curing (α ≈ 0%), a large number of ionized charge carriers in an alternating electrical field lead to the maximal relative permittivity Polymers 2021, 13, 3085 12 of 16 in the material. The number of ions reduces, and the ions’ mobility is restricted strongly by the increasing crosslink density, reflecting the reduced relative permittivity in a DEA measurement. The relative permittivity converges when the material almost finishes curing (α ≈ 100%), even though the small content of reaction agents that cannot find a reaction partner is still available. The regression results using Trust-Region algorithms for determining the unknown parameters regarding the DEA measurements are listed in Table 5. Polymers 2021, 13, x FOR PEER REVIEW 12 of 17 (a) (b) (c) (d) Figure 7. Reaction rate depending on crosslink density computed using varied cutoff distances of (a) 1.6 Å, (b) 1.7 Å, (c) 1.8 Å and (d) 1.9 Å. 3.3. Comparison between Simulation and Experimental Results 3.3.1. Reaction Kinetics Using DEA DEA measurements were performed under the pressure of 423 atm and at tempera- tures of 150 °C, 170 °C, and 190 °C. The change in reaction rate depending on the cross- linking amount, which is calculated according to the relative permittivity’s changes at dif- ferent temperatures, is plotted in Figure 8. At the beginning of curing (α ≈ 0%), a large number of ionized charge carriers in an alternating electrical field lead to the maximal relative permittivity in the material. The number of ions reduces, and the ions’ mobility is restricted strongly by the increasing crosslink density, reflecting the reduced relative per- mittivity in a DEA measurement. The relative permittivity converges when the material almost finishes curing (α ≈ 100%), even though the small content of reaction agents that cannot find a reaction partner is still available. The regression results using Trust-Region algorithms for determining the unknown parameters regarding the DEA measurements are listed in Table 5. Figure 7. Reaction rate depending on crosslink density computed using varied cutoff distances of (a) 1.6 Å, (b) 1.7 Å, (c) 1.8 Å and (d) 1.9 Å. Table 5. Parameters in the reaction model estimated by regression with a confidence interval of 95 %. Method Expression of Reaction Model f(α) cα nα Kα logA [1/s] R2 RMSE MD nth-order catalytic (nOk): f (α) = cα· (1− α)nα ·(1 + Kα·α) 0.3218 3.392 28.03 7.57 0.99 0.03 DSC 0.03 1.25 150.4 7.44 0.99 0.05 DEA 0.03 1.41 145.8 7.45 0.97 0.03 Polymers 2021, 13, 3085 13 of 16 Polymers 2021, 13, x FOR PEER REVIEW 13 of 17 Figure 8. Comparison of the kinetic model derived using MD, DEA, and DSC methods. Table 5. Parameters in the reaction model estimated by regression with a confidence interval of 95 %. Method Expression of Reaction Model f(α) 𝒄𝜶 𝒏𝜶 𝑲𝜶 logA [1/s] 𝑹𝟐 𝑹𝑴𝑺𝑬 MD nth-order catalytic (nOk): 𝑓(𝛼) = 𝑐 ∙ (1 − 𝛼) ∙ (1 + 𝐾 ∙ 𝛼) 0.3218 3.392 28.03 7.57 0.99 0.03 DSC 0.03 1.25 150.4 7.44 0.99 0.05 DEA 0.03 1.41 145.8 7.45 0.97 0.03 3.3.2. Reaction Kinetics Using DSC Concerning DSC measurements, the unknown parameters in Equation (8) are deter- mined in two steps. First, the reaction model f(α) is fitted with the Trust-Region algorithm, whereby the results are listed in Table 5. Figure 8 shows the fitted curve regarding the nth- order catalytic reaction model to the reaction rate, resulting in the dependency of the crosslinking amount measured with different heating rates β. Second, the Kissinger method was applied for determining the activation energy EA and the pre-exponential term A. Equation (9) can be evaluated according to a linear relationship 𝑦 = 𝑎 ∙ 𝑥 + 𝑏 , reported by Yan et al. [38]. Then, the activation energy EA and the pre-exponential term A are determined using Equations (10) and (11): 𝑙𝑛 𝛽𝑇 = − 𝐸𝑅 ∙ 𝑇 + 𝑙𝑛 − 𝐴 ∙ 𝑅𝐸 ∙ 𝑓 (𝛼 ) (9) 𝑎 = − => 𝐸 = 77.66 𝑘𝐽/𝑚𝑜𝑙 (10) 𝑏 = 𝑙𝑛 − ∙ ∙ 𝑓′(𝛼 ) => 𝐴 = ∙∙ 1/𝑠 (11) where 𝑓 (𝛼) = 𝑑𝑓(𝛼) 𝑑𝛼⁄ , Rgas is the gas constant, β is the heating rate, Tp represents the temperature, and αp is the crosslink density with the index p, which means that the indi- cated values correspond to the position of the rate peak maximum. Figure 8. Comparison of the kinetic model derived using MD, DEA, and DSC methods. 3.3.2. Reaction Kinetics Using DSC Concerning DSC measurements, the unknown parameters in Equation (8) are deter- mined in two steps. First, the reaction model f(α) is fitted with the Trust-Region algorithm, whereby the results are listed in Table 5. Figure 8 shows the fitted curve regarding the nth-order catalytic reaction model to the reaction rate, resulting in the dependency of the crosslinking amount measured with different heating rates β. Second, the Kissinger method was applied for determining the activation energy EA and the pre-exponential term A. Equation (9) can be evaluated according to a linear relationship y = aK·x + bK, reported by Yan et al. [38]. Then, the activation energy EA and the pre-exponential term A are determined using Equations (10) and (11): ln ( β TP2 ) = − EA Rgas·TP + ln ( − A·Rgas EA · f ′(αP) ) (9) aK = − EA Rgas => EA = 77.66 kJ/mol (10) bK = ln ( − A·Rgas EA · f ′(αP) ) = > A = −10bK ·EA Rgas· f ′(αP) 1/s (11) where f ′(α) = d f (α)/dα, Rgas is the gas constant, β is the heating rate, Tp represents the temperature, and αp is the crosslink density with the index p, which means that the indicated values correspond to the position of the rate peak maximum. 3.3.3. Comparison The results of the reaction model determined by different methods (MD, DSC, and DEA) are transmitted together in the Figure 8 (normalized dα/dt vs. α) and Table 5. While the maximum response rate appears at a low network level (approx. 20%) in the MD simulation, the maximum reaction rate occurs at a crosslink density of 50%, concerning ex- perimental results. According to the MD calculation, the crosslinking reaction ends almost at a high crosslink density of >85%. The mobility of free monomers/molecules is high at a low crosslink density, leading to the increased probability of the crosslinking reaction occurring. The barrier that prevents the molecular movements through morphological structures by the growing molecules is proportional to the increasing crosslink density. In an experimental investigation like DSC, the reaction heat is measured to determine the Polymers 2021, 13, 3085 14 of 16 crosslink density. The maximal crosslink density determined using experimental methods is 100% since it assumes that all monomers/functional groups meet a reaction partner. Additionally, the phase transition from solid to liquid during the melting phase is not considered to be using molecular dynamics. Therefore, after the phase change in the material, the movement of the reaction partners increases in the fluid state. Furthermore, the mixture of pure resin and hardener is modeled in molecular dynamics, without con- sidering the additives’ molecules which are commonly applied in a polymer blend for influencing the reaction speed or material properties. However, the comparison shows a good comparability and match between simulation and experimental results. 4. Conclusions In this work, molecular dynamics simulation was applied to depict the hardening process of an epoxy resin at the atomic level using a fixed force field PCFF under periodic boundary conditions with constant pressure and at different temperatures. The reactivity of the functional groups is varied through changing cutoff distance. In the first step, the influence of cutoff distance (1.48 Å→2.5 Å) and temperature (150 ◦C→190 ◦C) on the material’s curing is investigated. The numerical results show that the maximal achieved crosslink density increases with the growing cutoff distance during simulation time. The crosslink density increases up to a maximum value (e.g., approx. 90% at a cutoff distance of 1.8 Å) and then remains almost constant, demonstrating that a fully crosslinked state can hardly be reached. Additionally, the enhanced molecules’ diffusions due to increasing the temperature from 150 ◦C to 190 ◦C is significantly recognized in the curing computed with a low cutoff distance (e.g., 1.6 Å). However, the influence of rising environment temperature on the crosslink density is not significantly identified regarding numerical results computed with a cutoff distance of ≥ 1.7 Å. This is because the molecules have a higher probability of reacting with the neighboring reaction partners in a homogenized polymer mixture and the curing time is too short for the diffusion process, respectively. In the further step, the reaction model f(α), as a part of reaction kinetics, was determined based on the numerical results of the molecular dynamics study and compared to the experimental results determined by the DSC and DEA method. This approach indicates that the difference in the time scale between molecular dynamics and experiment need not be considered by adjusting the reaction cutoff distance and accepting accelerated reaction dynamics. The comparison demonstrates a similar course of the reaction model between simulation and experiment, while the maximum reaction rate in the MD simulation shifts to a lower crosslink density than the experiment. Nevertheless, the results show the feasibility of using molecular dynamics to derive reaction kinetics for a simple polymer system. In future investigations, the influences of chosen cutoff distances on the topology of the generated epoxy networks, and therefore the associated material’s mechanical properties, should be analyzed in detail. Aside from predicting the material properties, molecular dynamics simulation offers a possibility for estimating the maximum achievable crosslink density, which is expensive to measure in experiments. In addition, the determined reaction model combined with the activation energy can be implemented in a macroscale process simulation, e.g., using the FE method. Author Contributions: Conceptualization, S.S., H.Z. and S.Y.; methodology, S.Y. and W.V.; software, S.Y. and W.V.; validation, S.Y.; formal analysis, S.Y. and W.V.; investigation, S.Y. and W.V.; writing— original draft preparation, S.Y.; writing—review and editing, S.Y., W.V., H.Z. and S.S.; visualization, S.Y.; supervision, S.S. and H.Z. All authors have read and agreed to the published version of the manuscript. Funding: This research received no external funding. Institutional Review Board Statement: Not applicable. Informed Consent Statement: Not applicable. Polymers 2021, 13, 3085 15 of 16 Data Availability Statement: The data presented in this study are available on request from the corresponding author. Acknowledgments: The authors acknowledge support by the state of Baden-Württemberg thro- ugh bwHPC. Conflicts of Interest: The authors declare no conflict of interest. Appendix A 1 Figure A1. Schematic representation of the crosslinking reaction between (a: epoxide group and primary amine group, b: epoxide group and secondary amine group, and c: epoxide group and the hydroxyl group). References 1. 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. [CrossRef] 2. Gao, L.; Zhang, Q.; Li, H.; Yu, S.; Zhong, W.; Sui, G.; Yang, X. Effect of epoxy monomer structure on the curing process and thermo-mechanical characteristics of tri-functional epoxy/amine systems: A methodology combining atomistic molecular simulation with experimental analyses. Polym. Chem. 2017, 8, 2016–2027. [CrossRef] 3. Heine, D.R.; Grest, G.S.; Lorenz, C.D.; Tsige, M.; Stevens, M.J. Atomistic Simulations of End-Linked Poly(dimethylsiloxane) Networks: Structure and Relaxation. Macromolecules 2004, 37, 3857–3864. [CrossRef] 4. Li, C.; Coons, E.; Strachan, A. Material property prediction of thermoset polymers by molecular dynamics simulations. Acta Mech. 2014, 225, 1187–1196. [CrossRef] 5. Li, C.; Medvedev, G.A.; Lee, E.-W.; Kim, J.; Caruthers, J.M.; Strachan, A. Molecular dynamics simulations and experimental studies of the thermomechanical response of an epoxy thermoset polymer. Polymer 2012, 53, 4222–4230. [CrossRef] 6. 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. [CrossRef] 7. 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. [CrossRef] 8. Odegard, G.M.; Jensen, B.D.; Gowtham, S.; Wu, J.; He, J.; Zhang, Z. Predicting mechanical response of crosslinked epoxy using ReaxFF. Chem. Phys. Lett. 2014, 591, 175–178. [CrossRef] 9. Okabe, T.; Oya, Y.; Tanabe, K.; Kikugawa, G.; Yoshioka, K. Molecular dynamics simulation of crosslinked epoxy resins: Curing and mechanical properties. Eur. Polym. J. 2016, 80, 78–88. [CrossRef] 10. Schichtel, J.J.; Chattopadhyay, A. Modeling thermoset polymers using an improved molecular dynamics crosslinking methodology. Comput. Mater. Sci. 2020, 174, 109469. [CrossRef] 11. 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. [CrossRef] 12. Wu, C.; Xu, W. Atomistic molecular modelling of crosslinked epoxy resin. Polymer 2006, 47, 6004–6009. [CrossRef] 13. Yang, S.; Qu, J. Computing thermomechanical properties of crosslinked epoxy by molecular dynamic simulations. Polymer 2012, 53, 4806–4817. [CrossRef] http://doi.org/10.1016/j.polymer.2011.03.052 http://doi.org/10.1039/C7PY00063D http://doi.org/10.1021/ma035760j http://doi.org/10.1007/s00707-013-1064-2 http://doi.org/10.1016/j.polymer.2012.07.026 http://doi.org/10.1021/ma201390k http://doi.org/10.1016/j.polymer.2015.06.038 http://doi.org/10.1016/j.cplett.2013.11.036 http://doi.org/10.1016/j.eurpolymj.2016.04.019 http://doi.org/10.1016/j.commatsci.2019.109469 http://doi.org/10.1021/ma801153e http://doi.org/10.1016/j.polymer.2006.06.025 http://doi.org/10.1016/j.polymer.2012.08.045 Polymers 2021, 13, 3085 16 of 16 14. Huo, R.; Zhang, Z.; Athir, N.; Fan, Y.; Liu, J.; Shi, L. Designing high thermal conductivity of cross-linked epoxy resin via molecular dynamics simulations. Phys. Chem. Chem. Phys. 2020, 22, 19735–19745. [CrossRef] 15. Moeini, M.; Barbaz Isfahani, R.; Saber-Samandari, S.; Aghdam, M.M. Molecular dynamics simulations of the effect of temperature and strain rate on mechanical properties of graphene–epoxy nanocomposites. Mol. Simul. 2020, 46, 476–486. [CrossRef] 16. Li, J.; Jumpei, S.; Waizumi, H.; Oya, Y.; Huang, Y.; Kishimoto, N.; Okabe, T. A multiscale model for the synthesis of thermosetting resins: From the addition reaction to cross-linked network formation. Chem. Phys. Lett. 2019, 720, 64–69. [CrossRef] 17. Unger, R.; Braun, U.; Fankhänel, J.; Daum, B.; Arash, B.; Rolfes, R. Molecular modelling of epoxy resin crosslinking experimentally validated by near-infrared spectroscopy. Comput. Mater. Sci. 2019, 161, 223–235. [CrossRef] 18. Li, C.; Choi, P.; Sundararajan, P.R. Simulation of chain folding in polyethylene: A comparison of united atom and explicit hydrogen atom models. Polymer 2010, 51, 2803–2808. [CrossRef] 19. 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. [CrossRef] 20. Li, C.; Strachan, A. Molecular simulations of crosslinking process of thermosetting polymers. Polymer 2010, 51, 6058–6070. [CrossRef] 21. Doherty, D.C.; Holmes, B.N.; Leung, P.; Ross, R.B. Polymerization molecular dynamics simulations. I. Cross-linked atomistic models for poly(methacrylate) networks. Comput. Theor. Polym. Sci. 1998, 8, 169–178. [CrossRef] 22. Yarovsky, I. Computer simulation of structure and properties of crosslinked polymers: Application to epoxy resins. Polymer 2002, 43, 963–969. [CrossRef] 23. Diao, Z.; Zhao, Y.; Chen, B.; Duan, C.; Song, S. ReaxFF reactive force field for molecular dynamics simulations of epoxy resin thermal decomposition with model compound. J. Anal. Appl. Pyrolysis 2013, 104, 618–624. [CrossRef] 24. Pieter, J. EMC: Enhanced Monte Carlo. Available online: http://montecarlo.sourceforge.net/emc/Welcome.html (accessed on 31 May 2019). 25. Plimpton, S.; Roy, P.; Stevens, M. Fast Parallel Algorithms for Short-Range Molecular Dynamics. J. Comput. Phys. 1995, 117, 1–19. [CrossRef] 26. Plimpton, S. Particle-Mesh Ewald and rRESPA for Parallel Molecular Dynamics Simulations. In Proceedings of the Eighth SIAM Conference on Parallel Processing for Scientific Computing, Philadelphia, PA, USA, 14–17 March 1997. 27. Sandia National Laboratories. LAMMPS Molecular Dynamics Simulator. Available online: https://www.lammps.org/ (accessed on 15 April 2020). 28. Sun, H.; Mumby, S.J.; Maple, J.R.; Hagler, A.T. An ab Initio CFF93 All-Atom Force Field for Polycarbonates. J. Am. Chem. Soc. 1994, 116, 2978–2987. [CrossRef] 29. Gissinger, J.R.; Jensen, B.D.; Wise, K.E. Modeling chemical reactions in classical molecular dynamics simulations. Polymer 2017, 128, 211–217. [CrossRef] [PubMed] 30. Gissinger, J.R.; Jensen, B.D.; Wise, K.E. REACTER: A Heuristic Method for Reactive Molecular Dynamics. Macromolecules 2020, 53, 9953–9961. [CrossRef] 31. Sun, H. COMPASS: An ab Initio Force-Field Optimized for Condensed-Phase Applications Overview with Details on Alkane and Benzene Compounds. J. Phys. Chem. B 1998, 102, 7338–7364. [CrossRef] 32. Huheey, J.E. Anorganische Chemie: Prinzipien von Struktur und Reaktivität; De Gruyter: Berlin, Germany; New York, NY, USA, 1988; ISBN 3-11-008163-6. 33. Cîrcu, V.; Mocanu, A.S.; Roşu, C.; Manaila-Maximean, D.; Dumitraşcu, F. Thermal behaviour and electro-optical properties of a series of liquid crystals based on palladium complexes with mixed ligands. J. Therm. Anal. Calorim. 2012, 107, 877–886. [CrossRef] 34. Ros, u, C.; Manaila-Maximean, D.; Paraskos, A.J. Thermally Stimulated Depolarization Currents and Optical Transmission Studies on A 3,4-Dicyanothiophene-Based Bent-Rod Liquid Crystal. Mod. Phys. Lett. B 2002, 16, 473–483. [CrossRef] 35. Kissinger, H.E. Reaction Kinetics in Differential Thermal Analysis. Anal. Chem. 1957, 29, 1702–1706. [CrossRef] 36. Karcher, M.; Chaloupka, A.; Henning, F.; Schmölzer, S.; Moukhina, E. Bestimmung der Aushärtekinetik eines 2-stufig aushär- tenden Epoxidharzes zur Herstellung von Hochleistungsfaserverbunden. Z. Kunstst. 2016, 12, 80–111. [CrossRef] 37. Hardis, R.; Jessop, J.L.; Peters, F.E.; Kessler, M.R. Cure kinetics characterization and monitoring of an epoxy resin using DSC, Raman spectroscopy, and DEA. Compos. Part A Appl. Sci. Manuf. 2013, 49, 100–108. [CrossRef] 38. Yan, S.; Zeizinger, H.; Merten, C.; Schmauder, S. In-situ investigation of dielectric properties and reaction kinetics of a glass-fiber- reinforced epoxy composite material using dielectric analysis. Polym. Eng. Sci. 2021. [CrossRef] http://doi.org/10.1039/D0CP02819C http://doi.org/10.1080/08927022.2020.1729983 http://doi.org/10.1016/j.cplett.2019.02.012 http://doi.org/10.1016/j.commatsci.2019.01.054 http://doi.org/10.1016/j.polymer.2010.04.049 http://doi.org/10.1016/j.polymer.2007.02.007 http://doi.org/10.1016/j.polymer.2010.10.033 http://doi.org/10.1016/S1089-3156(98)00030-0 http://doi.org/10.1016/S0032-3861(01)00634-6 http://doi.org/10.1016/j.jaap.2013.05.005 http://montecarlo.sourceforge.net/emc/Welcome.html http://doi.org/10.1006/jcph.1995.1039 https://www.lammps.org/ http://doi.org/10.1021/ja00086a030 http://doi.org/10.1016/j.polymer.2017.09.038 http://www.ncbi.nlm.nih.gov/pubmed/33149370 http://doi.org/10.1021/acs.macromol.0c02012 http://doi.org/10.1021/jp980939v http://doi.org/10.1007/s10973-011-1609-3 http://doi.org/10.1142/S021798490200397X http://doi.org/10.1021/ac60131a045 http://doi.org/10.3139/O999.03022016 http://doi.org/10.1016/j.compositesa.2013.01.021 http://doi.org/10.1002/pen.25691 Introduction Materials and Methods Material Molecular Modeling Dielectric Analysis Differential Scanning Calorimetry Results and Discussion Curing Behavior Reaction Kinetics Using Molecular Dynamics Simulation Comparison between Simulation and Experimental Results Reaction Kinetics Using DEA Reaction Kinetics Using DSC Comparison Conclusions References