Kaiwen
Li
ab,
Gota
Kikugawa
*b,
Yoshiaki
Kawagoe
c,
Yinbo
Zhao
d and
Tomonaga
Okabe
cef
aDepartment of Finemechanics, Graduate School of Engineering, Tohoku University, 2-1-1 Katahira, Aoba-ku, Sendai 980-8577, Japan
bInstitute of Fluid Science, Tohoku University, 2-1-1 Katahira, Aoba-ku, Sendai 980-8577, Japan. E-mail: kikugawa@tohoku.ac.jp
cDepartment of Aerospace Engineering, Graduate School of Engineering, Tohoku University, 6-6-01, Aramaki Aza Aoba, Aoba-ku, Sendai, Miyagi 980-8597, Japan
dSchool of Aerospace Engineering and Applied Mechanics, Tongji University, Shanghai, 200092, PR China
eDepartment of Materials Science and Engineering, University of Washington, BOX 352120, Seattle, WA 98195-1750, USA
fResearch Center for Structural Materials, Polymer Matrix Hybrid Composite Materials Group, National Institute for Materials Science, 1-2-1 Sengen, Tsukuba, Ibaraki 305-0047, Japan
First published on 14th May 2024
The limitations in previous dissipative particle dynamics (DPD) studies confined simulations to a narrow resin range. This study refines DPD parameter calculation methodology, extending its application to diverse polymer materials. Using a bottom-up approach with molecular dynamics (MD) simulations, we evaluated solubility parameters and bead number density governing nonbonded interactions via the Flory–Huggins parameter and covalent-bonded interactions. Two solubility parameter methods, Hildebrand and Krevelen–Hoftyzer, were compared for DPD simulations. The Hildebrand method, utilizing MD simulations, demonstrates higher consistency and broader applicability in determining solubility parameters for all DPD particles. The DPD/MD curing reaction process was examined in three epoxy systems: DGEBA/4,4′-DDS, DGEBA/MPDA and DGEBA/DETA. Calculations for the curing profile, gelation point, radial distribution function and branch ratio were performed. Compared to MD data for DGEBA/4,4′-DDS, the maximum deviation in secondary reactions between epoxy and amine groups according to DPD simulations with Krevelen–Hoftyzer was 14.8%, while with the Hildebrand method, it was 1.7%. The accuracy of the DPD curing reaction in reproducing the structural properties verifies its expanded application to general polymeric material simulations. The proposed curing DPD simulations, with a short run time and minimal computational resources, contributes to high-throughput screening for optimal resins and investigates mesoscopic inhomogeneous structures in large resin systems.
Since there are numerous types of epoxy resin, the experimental selection of matrix resin is expensive and time-consuming. For this reason, molecular-scale curing simulations are expected to help reproduce and elucidate the thermoset resin reaction process. Molecular dynamics (MD) simulations reproduce molecular motions based on classical Newtonian mechanics. This has been proven to be an efficient way to perform curing reaction processes in various epoxy systems.8–12 However, it also requires considerable computational resources and has high time costs due to the large number of atoms and complex force field calculations in the MD system. Molecular-scale curing simulations are required for high-throughput screening of resins to assist in experimental material exploration. Therefore, saving computational time is a crucial consideration. Furthermore, most of the mesoscale structure-forming phenomena that take place in polymeric materials influence the material properties, such as reaction-induced phase separation.13 However, such mesoscopic structures cannot be captured in MD simulations. Hence, to meet the demands of high-throughput screening and accurately replicate mesoscopic inhomogeneous structures, there is a pressing need for a curing simulation capable of faster calculations for larger systems.
The dissipative particle dynamics (DPD) simulation method has gained increasing popularity in polymer simulations in recent years. DPD treats multiple atoms as coarse-grained particles, referred to as beads. These beads interact with each other through various forces, typically including dissipative, random, and conservative forces. In DPD simulations, the internal degrees of freedom of the beads are neglected, and only the motion of the beads is considered. Extensive studies have demonstrated that the DPD simulation method is an efficient coarse-grained simulation technique capable of capturing longer time and length scales.14–17 To enhance the accuracy of DPD simulations in reproducing the curing reaction process, the calculation methods and selection of interaction parameters play a crucial role in the overall DPD model. This paper focuses on precisely addressing these aspects, aiming to improve the fidelity of DPD simulations for capturing the intricate curing dynamics of epoxy materials.
DPD simulation is a simulation method aimed at enhancing computational efficiency and system scalability by consolidating atomic groups into coarse-grained particles. In 1992, Hoogerbrugge and Koelman introduced DPD as a computational approach to studying fluid flow and hydrodynamics,18,19 with the theoretical foundation later established by Español and Warren.20 Groot and Warren21 developed the conventional parameterization scheme for DPD interaction, linking the Flory–Huggins theory to interaction parameters. Building upon the aforementioned DPD force calculation method, Lísal et al. proposed a reaction ensemble DPD method to explore polydisperse homopolymer systems and supramolecular di-block copolymer melt.22,23 While their focus was on replicating the reaction equilibrium of polymer systems rather than the polymerization process itself, their work illustrated the applicability of DPD for simulating polymerization. Liu et al. developed a curing reaction model using the same DPD method and applied it to surface-initiated polymerization.24 However, in realistic systems, there are more than two kinds of beads, and the atomic groups represented in beads will influence their properties; specifically, the densities of different types of beads and the “bead–bead” interactions between different kinds of beads are variable.
There are several important approaches to parameterize DPD interactions for systems with variable bead–bead repulsions. Travis et al.25 constructed bead–bead DPD interactions based on bead solubility parameters and defined pair–wise interactions using the Hildebrand regular solution theory. However, in their work, bead–bead interactions are computed independently without considering how to recover correct bead densities. Alternatively, Backer et al.26 and Spaeth et al.27 introduced beads with variable masses at different coarse graining levels. In their study, the system consisted of three regions: one comprising “normal” DPD particles, another comprising larger DPD particles, and a small region at the intersection where both types of particles coexist. Interactions between the beads, known as “bead–bead” interactions, vary depending on the composition of these two types of beads. However, the densities of the “normal” beads remain constant throughout the system. Kacar et al.15 introduced a density-dependent parameterization allowing the recovery of experimental bead densities from simulations. In their method, they kept most of the Groot–Warren framework intact, including the Flory–Huggins theory, but extended it by allowing bead-type dependent densities and like–like parameters. They combined this DPD simulation method with MD simulation in a realistic epoxy system, using MD simulation for parameterization.28 In Kacar's framework, the calculation of solubility parameters utilized the Krevelen–Hoftyzer method,29 relying on directly deriving molar contributions from functional groups. However, for certain functional groups, data in Krevelen's dataset is unavailable, as is the case with the sulfur-containing group in 4,4′-DDS. Therefore, there is a need to expand the methodology for calculating DPD parameters proposed by Kacar et al. and apply it to a wider range of polymer materials for research.
For the curing reaction model, Komarov et al.30 proposed a coarse-grained curing molecular dynamics simulation to replicate the crosslinking process of thermoset resins. Li and Strachan31 applied this approach in MD simulation and used a novel method to update the charge after forming the crosslinked structure. Numerous studies have simulated curing processes using a distance-based criterion and constant reaction probability.32,33 However, these methods are not very accurate in reproducing the reaction process, as the reaction probability should vary depending on the type of reaction substance and reaction conditions. Therefore, a new reaction model was proposed by Okabe et al., which was applied to MD simulation.8 There are two criteria for the reaction process: the distance between the reaction sites and the Arrhenius-type reaction probability determined from the activation energy and local temperature. Then, Kawagoe et al.34 performed the curing reaction process by combining Okabe's curing reaction model with DPD simulation and introduced a reverse mapping procedure to reconstruct an all-atom MD system from a DPD system to verify the validity of the proposed curing DPD simulation and to evaluate its structural and thermomechanical properties. However, they employed the DPD scheme only for the single epoxy system given by Kacar et al.15 without discussing the calculation method of the interaction parameters. Consequently, to extend the application of the DPD method to a broader range of epoxy systems, it becomes essential to undertake a bottom-up exploration of DPD parameter calculation methods.
In order to effectively capture mesoscale structures and reduce computational costs, the DPD method has gained significant traction in the field of crosslinked polymers. Building upon previous research, the primary objective of this study is to develop a bottom-up approach for calculating DPD parameters, enabling their broader application in extensive epoxy systems. The DPD parameters are derived through MD simulations, primarily providing solubility parameters and bead densities, which are crucial for nonbonded interactions. Additionally, parameters related to covalent-bonded interactions, such as equilibrium distance and equilibrium angle, are also determined. To validate the accuracy of the developed DPD parametrization, structural properties are calculated in the on-reaction and post crosslinking process and compared against results obtained from MD simulations. By ensuring the fidelity of the DPD parametrization, this research aims to facilitate the widespread utilization of the DPD method in investigating complex epoxy systems and accurately reproducing their structural properties.
![]() | (1) |
![]() | (2) |
![]() | (3) |
![]() | (4) |
![]() | (5) |
![]() | (6) |
![]() | (7) |
rDPD3 = 3ρ−1. | (8) |
In this equation, δ can be evaluated directly from Hildebrand's theory,35 referred to as the Hildebrand solubility parameter hereafter. In the present study, we also employed the Krevelen–Hoftyzer's29 method, which offers an easier calculation based on the group contribution method to determine the solubility parameters. The two calculation methods for the solubility parameter mentioned above will be illustrated in the next section.
The remaining two forces FDij and FRij are given by
FDij = −γwD(rij)(![]() ![]() | (9) |
FRij = σwR(rij)θij![]() | (10) |
Español and Warren20 demonstrated that the system samples the canonical ensemble and adheres to the fluctuation–dissipation theorem if the following relationships hold:
wD(rij) = [wR(rij)]2, | (11) |
σ2 = 2γkT, | (12) |
![]() | (13) |
Travis et al.25 constructed bead–bead DPD interactions using bead solubility parameters and defined pair-wise interactions based on Hildebrand's regular solution theory. In Hildebrand's theories, the solubility parameter can be obtained using the following equations:
![]() | (14) |
Ecoh = ΔHvap − kT, | (15) |
ΔHvap = −Enonb + kT. | (16) |
![]() | (17) |
The above describes the calculation method for the Hildebrand solubility parameter. In numerous DPD studies, the Krevelen–Hoftyzer solubility parameter is also commonly employed to calculate the conservative force.15 In the estimation of the solubility parameter, the cohesive energy is decomposed into three components: the contribution of dispersion forces, Ed, the contribution of polar forces, Ep, and the contribution of hydrogen bonding, Eh. The cohesive energy is calculated as follows:
Ecoh = Ed + Ep + Eh. | (18) |
Accordingly, the calculation formula for the solubility parameter is as follows:
δ2 = δd2 + δp2 + δh2, | (19) |
The available experimental data, however, prove that it is impossible to derive a simple system for an accurate prediction of solubility parameter components from the chemical structure. In particular, the interaction of different structural groups in producing overall polar and hydrogen-bonding properties is so complicated that the solubility parameter components do not obey simple rules. Therefore, the approach of Hoftyzer and Van Krevelen29 was used in this study as follows:
![]() | (20) |
![]() | (21) |
![]() | (22) |
Structural group | F di (MJ m−3)0.5 mol−1 | F pi (MJ m−3)0.5 mol−1 | E hi (J mol−1) |
---|---|---|---|
–CH3 | 420 | 0 | 0 |
–CH2– | 270 | 0 | 0 |
![]() |
80 | 0 | 0 |
–NH2 | 280 | 0 | 8400 |
–NH– | 160 | 210 | 3100 |
![]() |
20 | 800 | 5000 |
–OH | 210 | 500 | 20![]() |
–O– | 100 | 400 | 3000 |
–C6H5 | 1430 | 110 | 0 |
![]() | (23) |
The friction coefficient γ of dissipative force is related to noise amplitude of random force as shown in eqn (23). Therefore, we only need to select an appropriate noise amplitude σ and time step Δt to determine the magnitude of dissipative force and random force.
The conversion from dimensionless to dimensional quantities in temperature in the DPD simulation utilizes the following formula:
![]() | (24) |
The time step Δt in the simulations is selected as 0.02tDPD. For our purpose the exact mapping of the time scale is not important since our main interest is the final crosslinked structure. Therefore, the unit DPD time, tDPD, can be calculated using the following dimensional calculation.
![]() | (25) |
We chose the DPD length scale such that the overall dimensionless DPD density is approximately ρrDPD3 = 3. This is the most used value in DPD simulations.21 To compute rDPD, the ‘real’ density ρ of the system needs to be known. This value can be determined from the experimental value of the number density. Here, we estimate the average molecular volume by assuming that this volume is the weighted average of the pure species volumes. The volume is calculated as follows:
![]() | (26) |
Referring to Groot,21 the dimensionless value of noise amplitude σ* was set to 3. The noise amplitude used in this study can be obtained by the following dimensional calculation.
![]() | (27) |
After determining rDPD, we can calculate the value of noise amplitude σ in the three systems which will greatly influence the dissipative force and random force.
For intramolecular interactions, FBij and FAijk are the two-body interaction of a bond and the three-body interaction of an angle, respectively.
VBij = kBij(rij − r0)2, | (28) |
VAijk = kAijk(θijk − θ0)2, | (29) |
Kacar et al.15 explored the complete construction method of DPD interaction parameters, including the discussion of the above bond and angle force constants. They found that all bead sequences form highly stiff bonds and stiff angles with different spring constants for different bond and angle species. However, these values were impractically high to be used in the DPD simulations. Therefore, relatively larger force constants were employed for all bond and angle species in this study, as follows, referring to Kacar et al.15
kBij = 500(kT/rDPD2), | (30) |
kAijk = 50kT. | (31) |
In this study, DPD parameters are obtained using MD simulation. This section introduces the simulation processes for the solubility parameter and number density of different types of beads, contributing to conservative force. Furthermore, the simulation process for the equilibrium distance and angle, which contribute to bonded force, is also described here.
First, the Winmostar36 software was utilized to visually construct both the full atomistic models corresponding to DPD beads and bead sequences. Subsequently, the Gaussian 16 C.0137 software was employed for the optimization of the structure (using the OPTIMIZE method) to obtain the molecular structure with the minimum energy. After breaking up an original molecule into fragments, the charges of the outermost atoms in the beads are retained as the charges of the original molecules. To ensure a net charge of 0, any excess charge is evenly distributed among the inner atoms, primarily the carbon atoms. To maintain the accuracy of simulating chemical reactions, the charges of oxygen and nitrogen atoms in the beads involved in the reaction groups are kept consistent with the charges of the original molecules. The assigned charges to each atom in the beads are shown in Fig. 2.
Due to the investigation of a large system and a shorter calculation time, the DREIDING force field38 was adopted. In the following procedures, we utilized the Nosé–Hoover chains thermostat39 to control the temperature and employed the Martyna–Tobias–Klein method40 in the NPT ensemble. The cut-off distance of 12 Å was adopted. Periodic boundary conditions (PBCs) were applied in all three dimensions. The r-RESPA multiple-time-step method41 was employed, with 1 fs for intermolecular forces and 0.2 fs for intramolecular forces. The smooth particle mesh Ewald (SPME) method42 was utilized to calculate the long-range Coulomb force in the in-house code.
To calculate solubility parameters and number density, firstly, the initial simulated box size was set at 50.0 × 50.0 × 50.0 Å3, with approximately 5000 atoms depending on the species of atomic groups. The temperature was set at 300 K, and the pressure at 1 atm, creating a liquid phase simulation system. Notably, C, H′′, I′, and F particles exist in the gaseous state under normal pressure, and to transition them to the liquid state, a minimum pressure was applied. Subsequently, 2000 step energy minimization using the steepest descent method was carried out. Following the energy minimization, a 1 ns NVT (constant number of particles, volume, and temperature) ensemble simulation and a 2 ns NPT (constant number of particles, pressure, and temperature) ensemble simulation were performed. The volume and nonbonded intermolecular energy could be directly obtained in the MD simulation.
To calculate the equilibrium distance and angle, we initially set the simulated initial box size at 230.0 × 230.0 × 230.0 Å3, with 30 bead sequences. The temperature was set to 300 K, and the pressure was maintained at 1 atm, establishing the simulation system as a gas phase. Subsequently, we initiated the simulation process with 2000 step energy minimization using the steepest descent method. Following this, a 2 ns NVT ensemble simulation was conducted. The equilibrium distance and angle were then determined by calculating the center of mass of the beads. By obtaining the probability distribution of bonds and angles, we traced the peak values, representing the most probable values, to ascertain the equilibrium distance and angle. These equilibrium values are utilized for initial state setup and in calculations of angle and bond potentials as shown in the formulas of eqn (28) and (29).
Epoxy system | Number of molecules | Average volumes (Hildebrand/Krevelen) (Å3) | r DPD (Hildebrand/Krevelen) (Å) | Noise level σ (Hildebrand/Krevelen) (10−17/J s0.5 m−1) |
---|---|---|---|---|
DGEBA/4,4′-DDS | 600/300 | 167.2/126.9 | 7.951/7.243 | 4.547/4.760 |
DGEBA/MPDA | 600/240 | 152.2/123.6 | 7.702/7.178 | 4.618/4.781 |
DGEBA/DETA | 600/300 | 164.9/120.5 | 7.902/7.121 | 4.557/4.801 |
To obtain the initial non-crosslinked structure, a relaxation step is necessary before the curing reaction. In the DPD simulation system, a DPD simulation of 0.693 ns which is equal to 10000 DPD simulation steps at 300 K was conducted prior to crosslinking. In the MD system, a series of simulations, including 2 ns of NPT (300 K, 1 atm) and 0.5 ns of NVT (300 K), were performed before the crosslinking process. Both DPD and MD simulations in this study underwent sufficient mixing and achieved equilibrium states before crosslinking. During the crosslinking process, the distances between all epoxy groups and amino groups in the system need to be calculated. Based on these distances, a determination is made regarding whether a chemical reaction occurs between all pairs of functional groups within a threshold value, denoted as Rc. The value of Rc is set to 5.64 Å, four times the equilibrium C–N bond length,31 in both the MD and DPD curing reaction processes.34 The settings for equilibrium angles and equilibrium distances used in this study are based on calculations from an all-atom MD model. Therefore, the corresponding value of Rc remains consistent. In the reaction determination, the temperature around the functional group should be calculated first. Fig. 3 is a detailed flow chart of the curing reaction process. In this procedure, the reaction probability was calculated based on the following Arrhenius equation.43
![]() | (32) |
To determine the reaction probability in the reaction calculation, the Monte Carlo method45 was used here. In the reaction judgement, the reaction probability is compared with a random number in the range of 0 to 1. If the reaction probability is greater than the number, a new bond is formed.
For the MD system, the following is the equation for calculating the kinetic energy after reaction.
Kafter = Kbefore + Hf, | (33) |
In the DPD system, an adjustment factor β was added such that the localized temperature increment in a coarse-grained system due to Hf is almost identical to that of an MD system. The following is the equation for calculating the kinetic energy after the reaction.
Kafter = Kbefore + βHf, | (34) |
The activation energy and the heat of formation determined using formulas from ref. 32–34 have been referenced from data in the previous publication.44 These values were the activation energy and heat of formation precisely calculated using the global reaction route mapping (GRRM) method46 on the basis of ab initio quantum chemical (QC) calculations.
Adding Hf induces chain reactions as T increases, i.e., k increases in eqn (32). During the relaxation phase of the simulation, multiple-step relaxation cycles were conducted to relax and equilibrate the new topology after cross-linking, and preparing for the next cross-linking step. For the MD curing reaction relaxation, the DREIDING force field was selected. Structural relaxation was carried out using the steepest descent method for 1000 steps. Subsequently, a 15 ps NPT (453 K, 1 atm) ensemble simulation and a 1 ps NVT (453 K) ensemble simulation were performed. In the relaxation process following DPD crosslinking, a 2.5 ps DPD simulation was executed. The crosslinking structure was formed using an external curing code written in the Python language, while the relaxation calculations were conducted using the in-house DPD/MD package. Following all cross-linking cycles, we conduct a relaxation process lasting 10 ns NPT (300 K, 1 atm) ensemble to obtain a stable structure.
Tables 3 and 4 present the number density of various bead types. The number density of newly generated beads is smaller than that of the corresponding original beads, indicating that the volume of the newly generated beads is larger. For example, the volumes of beads D′ and D′′ are larger than that of bead D. This expansion in volume is attributed to the crosslinking process, where carbon atoms are added, leading to a reduction in hydrogen atoms in beads D, F, and H. The relatively larger volume of carbon atoms contributed to the overall volume increase. The calculation of the Hildebrand solubility parameter involves nonbonded energy and system volume. Conversely, for the Krevelen–Hoftyzer solubility parameter, the volume was sourced from the SciFinder-n database,47 which is smaller compared to the volume obtained from MD simulations for the Hildebrand solubility parameter. This may be attributed to the Dreiding force field, which often leads to a larger volume due to the selection of LJ parameters. The LJ parameters in DREIDING were mainly derived from Buckingham potential (exponential-6, X6) parameters but without direct derivation from physical properties like density. The three components of the molar attraction function were determined by referencing the table provided by Hoftyzer and Van Krevelen.29 This study explores the differences between Krevelen–Hoftyzer and Hildebrand solubility parameters.
Bead type (DGEBA) | Number density (Å−3) | Bead type (DDS) | Number density (Å−3) | Bead type (MPDA) | Number density (Å−3) | Bead type (DETA) | Number density (Å−3) |
---|---|---|---|---|---|---|---|
A | 0.0061 | D | 0.0057 | F | 0.0175 | H | 0.0078 |
B | 0.0058 | E | 0.0126 | G | 0.0065 | I | 0.0065 |
C | 0.0056 | D′ | 0.0048 | F′ | 0.0078 | H′ | 0.0065 |
A′ | 0.0067 | D′′ | 0.0044 | F′′ | 0.0065 | H′′ | 0.0063 |
I′ | 0.0063 |
Bead type (DGEBA) | Number density (Å−3) | Bead type (DDS) | Number density (Å−3) | Bead type (MPDA) | Number density (Å−3) | Bead type (DETA) | Number density (Å−3) |
---|---|---|---|---|---|---|---|
A | 0.0088 | D | 0.0066 | F | 0.0201 | H | 0.0128 |
B | 0.0069 | E | 0.0139 | G | 0.0073 | I | 0.0085 |
C | 0.0083 | D′ | 0.0055 | F′ | 0.0128 | H′ | 0.0085 |
A′ | 0.0092 | D′′ | 0.0048 | F′′ | 0.0085 | H′′ | 0.0071 |
I′ | 0.0071 |
The evaluated solubility parameters are presented in Table 5. It is observed that the Hildebrand solubility parameters in this study are consistently smaller than the Krevelen–Hoftyzer solubility parameters. This deviation is partly attributed to the difference in volume of beads between the two methods. The volume calculated using the Hildebrand method is generally larger than that calculated using the Krevelen–Hoftyzer method, as indicated in Tables 3 and 4, resulting in a smaller calculated solubility parameter. Based on the DPD results compared with MD simulations, the Hildebrand method demonstrates a higher level of consistency. Furthermore, due to the limited data available for the Krevelen method, it cannot comprehensively calculate the solubility parameters for all DPD particles. In contrast, the Hildebrand method theoretically has the potential to calculate the solubility parameters for all DPD particles, making it more promising for widespread use. It is worth emphasizing this point in Tables 3–5: H and F′, H′ and F′′, and I, H′′, and I′ are the same type of particles with identical compositions, only differing in naming convention. Therefore, the volumes and solubility parameters of these particles are the same.
Bead type (DGEBA, DETA) | Hildebrand solubility parameter ((J cm3)0.5) | Krevelen–Hoftyzer solubility parameter ((J cm3)0.5) | Bead type (4,4′-DDS, MPDA) | Hildebrand solubility parameter ((J cm3)0.5) | Krevelen–Hoftyzer solubility parameter15 ((J cm3)0.5) |
---|---|---|---|---|---|
A | 15.34 | 16.44 | H | 18.61 | 19.99 |
B | 18.58 | 24.54 | I | 14.37 | 15.87 |
C | 9.37 | 12.35 | H′ | 14.37 | 15.87 |
A′ | 13.19 | 23.54 | H′′ | 18.38 | 19.38 |
D | 21.01 | 18.51 | I′ | 18.38 | 19.38 |
E | 25.35 | (25.35) | F | 10.97 | (10.97) |
D′ | 18.96 | 18.22 | G | 22.34 | (22.34) |
D′′ | 18.78 | 19.95 | F′ | 18.61 | 19.99 |
F′′ | 14.37 | 15.87 |
The equilibrium distance is a parameter used to calculate the bond force between two beads, while the equilibrium angle is used to calculate the potential for angle bending between three atoms. In this study, the equilibrium distances and angles between different bead sequences were obtained from MD simulations, as presented in Tables 6 and 7, respectively. The equilibrium distance and angle curves are provided in the ESI.† A total of 300 sample data points were obtained from the MD simulation, and the equilibrium distance was determined as the peak value of the maximum occurrence probability. These equilibrium values are utilized for the initial state setup and in calculations of angle and bond potentials as shown in the formulas of eqn (28) and (29).
Bead sequence (DGEBA) | Equilibrium distance (Å) | Bead sequence (curing agent) | Equilibrium distance (Å) | Bead sequence (between DGEBA and curing agent) | Equilibrium distance (Å) |
---|---|---|---|---|---|
A–B | 4.713 | D–E | 3.725 | A–D | 4.841 |
B–C | 3.883 | F–G | 2.809 | A–F | 3.885 |
H–I | 3.302 | A–H | 4.247 |
Bead sequence (DGEBA) | Equilibrium angle (Å) | Bead sequence (curing agent) | Equilibrium angle (Å) |
---|---|---|---|
A–B–C | 112.4 | D–E–D | 153.2 |
B–C–B | 143.5 | F–G–F | 178.2 |
H–I–H | 154.8 |
For the total simulation duration, DPD simulations require less time, primarily because they involve fewer simulation steps compared to MD simulations and fewer particles in the system. The MD simulation ran for approximately 50 hours on 200 cores, whereas the DPD simulation ran for only about 20 minutes on 20 cores, which is roughly 1/150 of the MD simulation.
The calculation formula for the deviation value is |NDPD − NMD|/NMD where NDPD represents the number of reactions obtained by a DPD curing reaction, and NMD represents the number of reactions obtained by an MD curing reaction. For the DGEBA/4,4′-DDS system, the DPD curing reaction profiles with Hildebrand parameters are more consistent with the MD curing reaction curve. In the DPD curing process with the Hildebrand solubility parameter, there is excellent agreement with MD simulations, in terms of both primary and secondary reactions. However, in the DPD curing process with the Krevelen–Hoftyzer solubility parameter, the ratio of primary reaction beads consistently falls below that of the MD curing reaction and the DPD curing process with the Hildebrand solubility parameter. Conversely, the ratio of secondary reaction beads is consistently higher. At the maximum deviation point, occurring near a 40% curing rate for the DPD curing reaction (Hildebrand) and MD curing reaction, the deviation of the primary reaction is 2.2%, and the deviation of the secondary reaction is 1.7%. In contrast, for the DPD curing reaction with Krevelen–Hoftyzer parameters, the maximum deviation point occurs near a 50% curing rate, with a deviation of the primary reaction of 25.8% and a deviation of the secondary reaction of 14.8%.
For the DGEBA/MPDA system, both types of DPD curing reactions yield similar results, although they deviate significantly from the MD curing reaction. This discrepancy can be attributed to the notably high size ratio of DPD particles, which influences the accuracy of the DPD simulation. For instance, considering the reactive particle F in MPDA, it comprises only four atoms, leading to a substantial size difference between bead F and the B particle. In the DGEBA/DETA system, the profiles of DPD curing reactions with the Hildebrand solubility parameters closely resemble those of the MD curing reaction compared to the Krevelen–Hoftyzer approach. Considering all three systems, the results obtained from the Hildebrand method are closer to those of MD simulations. Hildebrand solubility parameters are parameterized based on energy and volume from MD simulations using the Dreiding force field, which is the same force field applied in the MD reaction simulations of this study. The correspondence between MD and DPD (Krevelen) methods is not good, because, in this case, solubility is parameterized based on literature values rather than MD simulations with the Dreiding force field.
To elucidate the impact of variations in the curing reaction process on physical properties, the gelation point was examined. Criteria for identifying the gelation point included the largest molecular weight (LMW) and the second largest molecular weight (SMW). In this study, a more analytical mass-based criterion, denoted as the reduced molecular weight (RMW), was employed. The RMW is defined as the molecular weight average of all the molecules in the system. It increases until the largest molecule starts to dominate, after which it declines.
The gelation point, corresponding to the maximum point of the RMW profile concerning the curing rate, is utilized to assess the gelation process. From Fig. 5, gelation points for different systems can be determined. For the DPD curing reaction with the Krevelen–Hoftyzer solubility parameter, the gelation points are as follows: DGEBA/4,4′-DDS: 60%, DGEBA/MPDA: 65%, DGEBA/DETA: 55%. In the case of the DPD curing reaction with the Hildebrand solubility parameter, the gelation points are as follows: DGEBA/4,4′-DDS: 55%, DGEBA/MPDA: 65%, DGEBA/DETA: 50%. For the MD curing reaction, the gelation points are observed as follows: DGEBA/4,4′-DDS: 55%, DGEBA/MPDA: 60%, DGEBA/DETA: 50%. In the DGEBA/4,4′-DDS system, the gelation point obtained from the DPD curing reaction with the Hildebrand solubility parameter aligns well with that obtained from the MD curing reaction. Furthermore, it also shows good agreement with the theoretical estimation by M. Livraghi et al., suggesting a gelation point of about 57 ± 6% based on percolation theory.48
![]() | ||
Fig. 5 RMW, SMW LMW curve in (a) DGEBA/4,4′-DDS, (b) DGEBA/MPDA, and (c) DGEBA/DETA. Blue lines represent LMW curves, green lines represent SMW curves, and red lines represent RMW curves. |
The ratio of branches at the molecular scale represents the internal molecular connectivity of the crosslinking network structure within a material. This ratio is determined by the number of monomers with each branching structure, divided by the total number of monomers. It has implications for thermophysical properties such as glass transition temperature and Young's modulus. In a crosslinked system, each basic resin can undergo a maximum of two reactions, and each curing agent can undergo a maximum of four reactions. Consequently, resin/curing agent molecules will form varying numbers of branches, including 0, 1, 2, 3, or 4 branches. The ratios of branches can be compared between DPD simulation and MD simulation. In the DPD system, the reaction times of A, D, F, and H beads in a basic resin/curing agent molecule are calculated, while in the MD system, the reaction times of the amino group and epoxy group in a basic resin/curing agent molecule are counted. The number of branches changes continuously during the curing reaction. Before the gel point of 55%, branch 1 and branch 0 dominate the overall crosslinking network. From Fig. 6(a), it can be observed that branch 1 is located at the end of the crosslinked polymer chain. Therefore, during the period from the onset of cross-linking to the gel point, the increase in the proportion of branch 1 divides the system into smaller blocks of polymer. After exceeding approximately 40% conversion, the proportion of branch 1 starts to decrease, indicating that the dispersed polymer blocks begin to coalesce into larger polymer structures. Until reaching the gel point of 55%, the largest weight polymers begin to form and continue to expand with increasing conversion.
In the DGEBA/4,4′-DDS system, the ratio of branches in the DPD curing process, using the Hildebrand solubility parameter, closely aligns with that observed in the MD curing reaction process. This similarity arises from the comparable ratios of primary and secondary reactions observed in the MD system, as compared to the DPD curing process using the Hildebrand solubility parameter. As a result, the branch structures remain relatively similar during the curing reaction process. However, in the DPD curing process with the Krevelen–Hoftyzer solubility parameter, the growth rate of primary reactions increases significantly after reaching a 20% curing rate. Consequently, during this stage, secondary reactions occur slowly following the primary reactions. For example, the generation of 1 branch from the primary reaction quickly occurs, but the change to 2 branches is slow. As a result, the curve representing the presence of 1 branch reaches a peak in the DPD curing process, whereas the curve for 1 branch in the MD curing process represents a plateau.
Afterwards, the DPD/MD curing reaction process was conducted in three epoxy systems: DGEBA/4,4′-DDS, DGEBA/MPDA and DGEBA/DETA. Subsequently, we compared the number of reactions for different types. We observed that the DPD simulation with the Hildebrand solubility parameter is more consistent with MD simulations compared to the DPD simulation with the Krevelen–Hoftyzer solubility parameter. The Curing properties in the reaction and the structure properties were compared between these methods. As a result, the DPD simulation can reproduce the structure of the full atomistic model to a certain degree. From a comprehensive perspective of the above structural parameters, the Hildebrand parameter is more suitable for reproducing the curing process in MD simulations. In comparison to the Krevelen–Hoftyzer method, the Hildebrand method, based on comprehensive simulations to calculate solubility parameters, theoretically enables the calculation of solubility parameters for any type of DPD particles. This conclusion significantly expands the applicability of DPD simulations for polymeric materials.
Footnote |
† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d3sm01743e |
This journal is © The Royal Society of Chemistry 2024 |