Characterization of Cure Behavior in Epoxy Using Molecular Dynamics Simulation Compared with Dielectric Analysis and DSC

The curing behavior of a thermosetting material that influences the properties of the material 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.


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 experimental 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 material design and to predict material properties. The previous studies [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15] focused on predicting the thermomechanical properties (e.g., glass transition temperature, Young's modulus, thermal conductivity, and compressibility) in already cured molecular structures of the epoxy polymers using MD simulation. The molecular structures (Li et al. [4]) 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 (R 1 -CCO) and the primary amine group (R 2 -NH 2 ); (b) Reaction between the epoxide group (R 1 -CCO) and the secondary amine group (R 3 -NH-R 4 ); (c) Reaction between the epoxide group (R 1 -CCO) and the hydroxyl group (R 5 -OH). The detailed sequence of the crosslinking reaction with molecular structures as examples, is shown in Figure A1 in Appendix A.

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][20][21][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).

•
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/cm³.
As an initial equilibration, a short simulation sequence in microcanonical ensemble 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.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 examples, is shown in Figure A1 in Appendix A.

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][20][21][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).

•
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/cm³.
As an initial equilibration, a short simulation sequence in microcanonical ensemble 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.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 examples, is shown in Figure A1 in Appendix A.

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][20][21][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).

•
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 N BADGE = 1817, which included the amount of the epoxide group N -CCO = 3634. The amount of the hardener EN N EN = 908, and includes the primary amine group N -NH2 = 1818. The initial density after creation by EMC was 1 g/cm 3 .
As an initial equilibration, a short simulation sequence in microcanonical ensemble 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 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.

•
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 D cut , 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 D cut . 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 D cut ), 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 V bond between the atom i and j were generated with Equation (1) using class2 bond style presented by Sun [31]: Parameter r ij.0 presents the equilibrium bond distance. K 2 , K 3 , and K 4 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 D cut on the cure kinetics were analyzed. The curing process was simulated at temperatures of 150 • C, 170 • C, and 190 • C while varying the cutoff distance D cut of the O-H bond from 1.5 Å→2.5 Å.

Dielectric Analysis
The dielectric analysis is a measuring procedure to study the material's response to an applied electric field. The electrical current I electr (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 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 I electr 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).
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 t 0 , 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 t 0 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 t 1 when the relative permittivity's change rate is less than 1%. Thereby, the value of relative permittivity ε r,1 is determined at time t 1 . 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 is the crosslink density determined using DEA. The term ε r,0 means the maximal relative permittivity occurring at t 0 ; ε r,1 is the relative permittivity determined at time t 1 , when the reaction is terminated, and ε r,∆ is the value of the difference between ε r,0 and ε r,1 .

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 N 2 as the purge gas.

Differential Scanning Calorimetry
DSC is a classic method for analyzing the correlation between the curing amou temperature, and time, based on the amount of heat released from crosslinking reactio The equipment DSC Q100 (TA Instruments-Waters LLC, New Castle, DE, USA) was u for studying the reaction kinetics. Non-isothermal measurements were performed w heating rates of 2 °C/min, 5 °C/min, 10 °C/min, 15 °C/min, and 20 °C/min with nitro N2 as the purge gas.
The thermally stimulated depolarization currents method [33,34] and DSC are ap cable for an experimental determination of the activation energy, while the DSC meth has been employed in our study. Concerning the results of DSC measurements, the a vation energy EA and the pre-exponential A could be determined using the Kissin method [35]. For predicting the unknown parameter in reaction model f(α), the follow procedure is performed. During the curing process, the change of heat flow in depe ency of time and temperature is detected by DSC shown in Figure 2a. According to resultant change of heat amount, temperature-and time-dependent crosslink density ratio of the released heat amount Hrelease to the maximum heat amount Hmax in a fully cro linked state, using Equation (4) (see Figure 2b). Figure 2c shows the reaction rate in pendency on the crosslinking degree to estimate the unknown parameter in the react model f(α): The thermally stimulated depolarization currents method [33,34] and DSC are applicable 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 E A 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 H release to the maximum heat amount H max 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(α): The term α DSC represents the crosslink density determined using DSC. The H release is the released heat amount by the material in the partially cured state, and H max is the maximum heat amount released by the material in the fully cured state.

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

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 N u,CCO to the total number of epoxide groups N total,CCO .
The computing time and effort, with a calculated bond dissociation length (D cut = 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 D cut 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, 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]. is set to larger than 1.6 Å, the crosslinking process becomes significantly faster. Additionally, 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]. 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.
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 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.
The term α MD represents the crosslink density determined using molecular dynamics simulation. N u,CCO is the number of epoxide groups already utilized during curing, and N total,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. 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.
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 increasing temperature of 150 °C→190 °C, which correlates to the thermal expansion of the molecular system influenced by a higher temperature. 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.  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. 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 experiment contains additives such as the reaction inhibitor, which delays and slows down the curing reaction.
Additionally, monomers are homogeneously distributed in the simulation volume, lead ing to a shorter balancing time of the molecular movements and a faster reaction. Furthe more, besides pure epoxy resin and hardener, the commercial EP material investigated i the experiment contains additives such as the reaction inhibitor, which delays and slow down the curing reaction.

Reaction Kinetics Using Molecular Dynamics Simulation
Reaction kinetics describe the rate of the chemical reaction, which considers tempe ature, pressure and concentration dependency: k(T), f(α), and h(p) (see Equation (6)).
Supposing that the pressure p is kept constant during the chemical reaction proces the pressure dependency can be ignored, and Equation (6) is simplified into Equations (7 and (8).
The temperature dependency, k(T) = A·exp(-EA/(Rgas·T)), generally known as the A rhenius equation, consists of the material-related activation energy EA and the pre-expo nential term A: The reaction model f(α) describes the reaction mechanisms in dependency on the pro portion of crosslinked molecules. There are already empirical equations available (such a nth-order catalytic, Avrami-Erofeev, Prout-Tompkins, and Sestak-Berggren reactio model [36,37]), whose unknown factors can be determined by regression methods. Sinc the empirical model (nth-order catalytic) best fits the numerical results for this invest 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 eva

Reaction Kinetics Using Molecular Dynamics Simulation
Reaction kinetics describe the rate of the chemical reaction, which considers temperature, 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).
The temperature dependency, k(T) = A·exp(-E A /(R gas ·T)), generally known as the Arrhenius equation, consists of the material-related activation energy E A and the preexponential term A: 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 simulation 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 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 temperatures 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 temperature 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). uation 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 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 temperatures 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 temperature 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).  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 increasing 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.  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 increasing 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.

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 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.

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 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    Concerning DSC measurements, the unknown parameters in Equation (8) are determined 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 nthorder 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): 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 indicated values correspond to the position of the rate peak maximum.

Reaction Kinetics Using DSC
Concerning DSC measurements, the unknown parameters in Equation (8) are determined 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 E A and the pre-exponential term A. Equation (9) can be evaluated according to a linear relationship y = a K ·x + b K , reported by Yan et al. [38]. Then, the activation energy E A and the pre-exponential term A are determined using Equations (10) and (11): where f (α) = d f (α)/dα, R gas is the gas constant, β is the heating rate, T p 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.

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 experimental 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 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 considering 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.

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.