Properties prediction and design of self-healing epoxy resin combining molecular dynamics simulation and back propagation neural network

It is difficult to synergistically improve the self-healing and mechanical properties of self-healing epoxy resin system based on Diels-Alder (DA) reversible covalent bond due to its different compositions. In this study, a method combining molecular dynamics (MD) simulation and back propagation (BP) neural network was employed to achieve a balance between self-healing property and mechanical property. In addition, self-healing and mechanical properties were predicted and the influence of different compositions of the system on these properties was investigated by the BP neural network. It was found that the DA bond formed by the furan group of furfuryl glycidyl ether (FGE) and the maleimide group of bismaleimide (BMI) was of great importance for self-healing property, and that diglycidyl ether of bisphenol-A (DGEBA) with a mass fraction of 25%–48% could improve the tensile property of the system. Finally, a self-healing epoxy resin system with superior comprehensive properties was designed through the proposed method under four indexes of Young’s modulus (E), self-healing efficiency ( η ), Ultimate Tensile Strength (UTS), and glass transition temperature (Tg).


Introduction
Epoxy resin has been playing an important role in coatings, adhesives, electronic and electrical materials, engineering plastics and composite materials due to its advantages such as strong adhesion, good corrosion resistance, excellent mechanical property, simple coating process and low cost [1][2][3]. However, cracks may arise inside or on the surface of epoxy resin during long-term use, which will degrade its material performance and reduce its service life, thus increasing the potential safety hazards in the process of material use [4]. In addition, the macromolecules of the crosslinked three-dimensional network in the epoxy resin are insoluble and infusible [5]. It is also difficult to be degraded and decomposed by microorganisms. Therefore, its relevant wastes are difficult to treat and even cannot be treated or recycled. Its internal micro-cracks are difficult to be healed by conventional methods such as welding, cementing and repairing [6]. And its structure can be destroyed if recycling methods such as particle, chemical and energy recycling is adopted. More importantly, its recycling rate is low, which does not alleviate environmental pollution [7,8]. Therefore, in order to fully recover epoxy resin without changing its three-dimensional crosslinked structure, it is of great theoretical value and practical significance to realize its self-healing function.
Intrinsic self-healing epoxy resin is a thermosetting polymer, which has excellent self-healing property due to reversible chemical reaction under external light or thermal stimulation. The intrinsic self-healing epoxy resin based on Diels-Alder (DA) reversible covalent bond is superior among self-healing polymers. DA reaction is Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. characterized by strong reversibility, high yield, mild and simple reaction conditions, etc, and can achieve multiple self-healing of materials, so it can be used for the synthesis, structural design and performance evaluation of epoxy resin system. Peterson et al [9] synthesized a DA bond-containing epoxy resin by using diene-containing phenyl glycidyl ether and dienophilic bismaleimide. He et al [10] fabricated epoxy resin by the DA reaction between diglycidyl furfurlamine and bismaleimide and found that it exhibited outstanding multiple damage/self-healing performance. Athanasios et al [11] studied the effects of healing agent concentration and curing period on the self-healing efficiency of carbon fiber reinforced composites based on DA reaction. However, the proposed self-healing epoxy resin systems based on DA reversible covalent bond does not solve the contradiction between self-healing and mechanical properties. The effect on self-healing and mechanical properties varies among proportion of each component. In order to achieve a balance between self-healing and mechanical properties, we proposed a method to predict properties and design the components of self-healing epoxy resin system so as to improve its mechanical property while ensuring its self-healing property.
The combination of molecular dynamics (MD) simulations and back propagation (BP) neural network is fit for designing resin. The performance of resin systems with different compositions can be automatically predicted by the BP neural network when the sufficient data obtained by MD simulation is used as its raw data. On the one hand, MD simulation can analyze the molecular topological characteristics of polymers from the atomic level and fundamentally find the reasons for differences in resin property [12][13][14], thereby greatly improving the efficiency of designing resin components. Yang et al [15] developed an interactive crosslinking relaxation method to construct a simulation unit and studied the dependence of elastic modulus, Poisson's ratio and coefficient of thermal expansion on crosslinking density and temperature. Jin et al [16] studied the effects of epoxy components on mechanical and thermal properties and designed a high-performance epoxy resin system by MD simulation. On the other hand, as a kind of artificial neural network in machine learning (ML), BP neural network is a multi-layer feedforward network trained by error back propagation algorithm, which can extract knowledge from training data and identify its internal rules. And then these internal rules can be applied to new data for further prediction [17][18][19][20]. Kokin et al [21] predicted the physical properties of thermosetting resin by using ML method. Hussein et al [22] utilized BP neural network method to explore the factors affecting the transverse buckling strength of corrugated composite panel. The combination of MD simulation and BP neural network does not require massive simulation calculations and heavy experimental workload [23,24].
On this basis, we firstly determined the consistency between experimental results and MD simulation results and verified the feasibility of MD in the self-healing epoxy resin system. Secondly, the glass transition, stretching and self-healing of resin systems with different compositions were simulated by MD, and the simulated data was imported into BP neural network. Finally, the performance of the system was predicted by BP neural network and experimentally verified. A self-healing epoxy resin system with high performance was designed by this method. Additionally, the effect of component on the performance of self-healing epoxy resin system was investigated, and the healing mechanism of its microstructure was studied by MD simulation.

Curing experiment of epoxy resin
In this study, epoxy resins were fabricated by furfuryl glycidyl ether (FGE) and diglycidyl ether of bisphenol-A (DGEBA), and Methylhexahydrophthalic anhydride (MeHHPA) and Bismaleimide (BMI) were used as curing agents.
The whole reaction mechanism in the bulk matrix is schematically represented in figure 1. The DA reaction occurred between furan and maleimide groups, and then DA adduct and DGEBA were respectively subjected to epoxy ring-opening reaction to obtain two kinds of resins A and B [25,26]. In the self-healing process, resin A experienced reverse DA (r-DA) reaction under certain conditions to generate BMI and ring-opened FGE (r-FGE), and then the DA reaction of these two substances occurred under different conditions to regenerate resin A [27].

Materials
Epichlorohydrin, tetrabutylammonium hydrogen sulfate, furfuryl alcohol, ethylacetate and sodium hydroxide supplied by Shanghai Aladdin Biochemical Technology Co., Ltd were used as reagents for synthesizing epoxy monomer FGE. The curing agents MeHHPA and BMI were purchased from Guangzhou Hexing Chemical Co., Ltd The epoxy resin DGEBA was supplied by Changzhou Lebang Composite Material Co., Ltd All the reagents were of analytical grades and used without any further modification.

Synthesis of FGE
Epichlorohydrin (90 ml) and tetrabutylammonium hydrogen sulfate (4 g) were charged into a three-necked flask equipped with a thermometer, an agitator and an inlet for dry nitrogen. Furfuryl alcohol (86 ml) was slowly added dropwise to the system through a constant pressure funnel at the temperature of below 25°C under the protection of nitrogen. After reaction at room temperature for 4 h, 50% NaOH solution (160 ml) was added, and the reaction was continued at room temperature for 2 h. Then the obtained product was extracted with ethylacetate to obtain the upper liquid phase. After the ethylacetate completely volatilized, the liquid phase was washed with water for 3 times. The lower oil phase was collected and the residual moisture was removed by vacuum distillation to obtain an epoxy resin monomer FGE containing a furan ring.

Preparation of DA-based self-healing epoxy resin
BMI was dissolved in FGE and DGEBA and stirred for 10 min at 90°C to obtain a homogeneous system. Then, MeHHPA was added, and the mixture was stirred at 80°C to get a clear homogeneous solution. The resultant homogeneous compound was degassed in a vacuum desiccator at 60°C for 30 min and then poured into a preheating mold at 60°C. The curing was conducted for 3 h at 70°C, 3 h at 100°C, 3 h at 120°C, 3 h at 150°C, and 24 h at 70°C [28].

Characterization and testing
A small amount of resultant homogeneous compound was taken out for differential scanning calorimetry (DSC) test. DSC (TA Q20 type) tests were carried out following ASTM D3418-15 test standard at the temperature range of 25°C to 200°C (Heating rate 10°C min −1 , nitrogen atmosphere). The heat flow j of the curing process was measured, and then the conversion degree α was calculated according to the measurement results. The relation between α and j is shown in the following formula [29]: where ΔH is the enthalpy change of the chemical reaction.
Dynamic thermomechanical analysis (DMA, TA Q800 type) was performed according to ASTM E1640-18 test standard with a scanning range of 30°C to 200°C (Heating rate 1°C min −1 , frequency 1 Hz). Three-point bending mode was adopted in the experiment. The system storage modulus was measured and recorded. Two tangents to the storage modulus curve below the transition temperature and at the inflection point which was approximately midway between the sigmoidal changes associated with the transition were constructed. The temperature at which these tangent lines intersected was reported as glass transition temperature (Tg).
Tensile tests were executed on a computer-controlled electronic universal testing machine (WDW-100 type) in accordance with ASTM D638-10 test standard (Stretching rate 2 mm min −1 ).
The crack healing process was characterized by field-emission scanning electron microscopy (FE-SEM, JEOL JSM 6700F type) at an emission current of 12 μA and an accelerating voltage of 5 kV.
The self-healing performance was characterized by tensile test. The sample with tensile cracks was healed at 120°C for 30 min and then at 70°C for 36 h. The last step, the healed sample was tested again according to the above tensile test operation. Self-healing efficiency (h) was calculated by the percentage of healing tensile strength to the original tensile strength.

Molecular dynamics simulation for data collection
In the MD simulation, MeHHPA was always able to accurately open the epoxy ring of FGE and DGEBA epoxy resins. We designed 25 groups of self-healing epoxy resin systems (table 1) with different compositions by changing the ratio of FGE, BMI, and DGEBA. In each composition, a sample number was assigned to each mixture ratio, which is expressed by mole ratio.
Each stage of MD simulation was realized by COMPASS force field [30,31]. The molecular models of epoxy monomers FGE and DGEBA, and curing agents MeHHPA and BMI were created according to its structure. The mixture of epoxy monomers and curing agents with the required stoichiometric ratio was filled into the simulation cell with three-dimensional periodic boundary according to the ratios in table 1. The conjugate gradient method was used to minimize the energy, and the configuration with the minimum potential energy was selected. This configuration was equilibrated in an isothermal and isochoric (NVT) MD simulation for 60 ps at 600 K. Subsequently, MD simulation was executed under an isothermal and isobaric (NPT) for 100 ps at atmospheric pressure to eliminate the internal stress of the system. The temperature and pressure were equilibrated with Andersen and Berendsen algorithms [32,33], respectively.  In the curing process, crosslinking reactions were simulated in a stepwise mode by distance-based standard. Bonds were formed between atoms that reacted within a truncation distance determined by the van der Waals radius of the reacting atoms. As the crosslinking density increased, it became more difficult to form new bonds. The truncation distance was gradually increased to accelerate the process. The truncation distance of O-C bond formed by crosslinking was set between 3.5 Å and 7 Å and that of DA bond between 4 Å and 8 Å. In order to obtain the self-healing epoxy resin crosslinked model (figure 2), the preset conversion limit in the simulation was 90% in all cases [34].The temperature and pressure were equilibrated with Andersen and Parrinello-Rahman algorithms [33], respectively.
A cooling simulation with a cooling rate of 10 K/60 ps at 500K to 300K was performed on the above crosslinked model during the glass transition. The temperature and pressure were equilibrated with Nose and Parrinello-Rahman algorithms [33], respectively. In order to plot a density-temperature curve, the system density was recorded every time the temperature dropped by 25 K until it reached 300 K. The change in the slope of a density-temperature curve was expressed by Tg. The stretch simulation was performed by expanding the length of the cubic unit along the stretching direction at each dynamic procedure. The ultimate strain was set to 2%, and the strain rate was 5×10 8 s −1 . The temperature and pressure were equilibrated with NHL and Souza-Martins algorithms [32,33], respectively. The stress-strain curve was obtained by a series of stress-strain points in the tensile process. Young's modulus (E) was calculated by linearly fitting the straight segment of the stressstrain curve, and Ultimate Tensile Strength (UTS) was defined as the maximum stress obtained from the stressstrain curve.
Self-healing process was simulated after the crosslinking and stretching procedure described earlier. As shown in figure 3, the simulation of the self-healing process mainly includes three steps, the disconnection of DA bonds, the rearrangement of small molecular structures, and the generation of DA bonds. To begin with an equilibration of the damaged model, one-time equilibration with NVT and then NPT dynamics were performed for 65 ps and 55 ps respectively. Afterwards a 70 ps NPT MD simulation was conducted at a temperature of 120°C and a pressure of 0.3 MPa, followed by one-time NPT dynamics for 150 ps at atmospheric pressure and 120°C. This multi-step relaxation program slowly opened the DA bonds, two C-C single bonds in the sixmembered ring were broken at the same time. Of the remaining four C-C bonds, three C-C single bonds became C-C double bonds in the meanwhile, and the other C-C double bond was transformed into C-C single bond at this time. Subsequently, a geometry optimization was performed for 200 ps to reduce energy. The minimization and stabilization of energy lead to the rearrangement of small molecular structures in the model until the structures was in equilibration state. And an NVT dynamics was performed for 50 ps at 600 K to obtain an optimized model. Finally, the reaction sites were defined, DA bonds were formed between reaction sites that within a truncation distance by continuously updating the reaction radius. The formation of the DA bonds in each truncation distance was performed under an NPT ensemble at a temperature of 70°C for 100 ps. From 4 Å to 8 Å, five steps were carried out in sequence, a total of 500 ps bonding simulation. As truncation distance gradually increased, the bonding process stopped when the truncation distance reached the maximum value or the crosslinking degree reached the conversion limit. The butadiene and ethylene structure gradually formed a six-membered ring through the set truncation distance. The C-C single bonds and double bonds in the original butadiene and ethylene structures were mutually converted. Therefore, the final healed model was established. According to the above stretch simulation procedure, the final model was simulated to obtain mechanical properties.

Machine learning techniques
The data on self-healing and mechanical properties were collected from MD simulation results. The property prediction and design of self-healing epoxy resin were realized using BP neural network. Using gradient descent method, BP neural network adjusted the weight and threshold to continuously reduce the errors between the actual and expected output value and minimized the amount of data required for training the network with the same accuracy.  The architecture of the BP neural network model is depicted in figure 4. The model consists of an input layer, three hidden layers, and an output layer. The thermal and mechanical properties of polymers directly affect the functional and processing performances of the material. The decline of these properties determines their service life and service quality. The self-healing property of the polymer reflects its healing and recycling value. The molecular network structure of thermosetting polymers plays a critical role in these material properties, and the molecular network is determined by the composition of the polymer system. Therefore, the content of FGE, BMI, and DGEBA was used as the variables for the input layer in the epoxy resin system, and Tg, E, UTS, and h as the results of the output layer. The compositions and their MD simulation results of the 25 designed systems in table 1 were defined as raw data set. The operation details of the BP neural network are listed in table 2.
Variables were transferred in the input layer, hidden layer and output layer of the BP neural network. Assuming that the input vector is X, and the output vector is Y, the weight matrix and deviation matrix from the input layer to the hidden layer are W 1 and V 1 , and those from the hidden layer to the output layer are W 2 and V 2 , if the neuron transfer functions of the hidden layer and output layer are C 1 and C 2 respectively, then the output of the BP neural network is: The whole operation process of ML is described in figure 5, including three steps of training, test and prediction. The raw data set was divided into a training data set (80%) and a test data set (20%). The training data set was randomized to avoid data sorting. In order to improve the generalization capacity of BP neural network, the data set was normalized in the range of 0-1, which is expressed as: where X is the measured value, X min and X max are the minimum and maximum values, and N is the normalized value. The mean square error (MSE) was used as the objective function to evaluate the training effect of the BP neural network model in this work. The MSE is expressed by the following equation: where p(n) and t(n) are the predicted and target values, and L is the number of training data samples. The optimizer was improved by using Adam algorithm, which mainly solved the weight increment of traditional neural network iteration by moving average iteration. The over-fitting problem was well solved when the initial iterative error was large and the training sample was small. In the back propagation of the Adam network, we first tried to update the error using the stochastic gradient descent optimization algorithm and improve the error back propagation using the Adam algorithm. The neural network initial weights were presented and then adjusted by the loss value between the predicted value and the target value using the optimizer. As iteration of network training increased, the weight values were gradually optimized in the correct direction, and the corresponding MSE loss was progressively decreased. When the number of iterations in the training cycle was sufficient, the number of iterations to minimize the loss function could be obtained. The weights corresponding to the number of iterations made the output value of the network closest to the target value, forming a well-trained neural network.
The neural network model after being trained was verified by the test data set and all the predicted results were obtained. The results of more than 1,000 systems were predicted by adjusting the proportions of FGE, BMI, and DGEBA.

Results and analysis
4.1. Self-healing process DSC analysis was performed on the DA-based self-healing epoxy resin systems to study the DA reaction and the reverse DA (r-DA) reaction. It can be seen from figure 6 that two exothermic peaks were observed at around 70°C and 120°C, of which the former was attributed to the formation of DA adducts while the latter was attributed to the formation of r-DA bonds. The emergence of these two characteristic peaks indicated that DA and r-DA bonds were successfully formed in the curing process, thus making self-healing possible. As shown in figure 7, the microscopic images of the original resin pattern were intact and crack-free. After the sample was damaged by external force, the cracks on the sample were completely visible, which could be regarded as a local fracture from a microscopic view. With the occurrence of r-DA and DA reaction, these cracks gradually disappeared. Self-healing microscopic images showed no cracks at 100 and 200 exposures, and only shallow traces at 500 exposures, indicating that the self-healing epoxy resin based on DA reaction had an excellent self-healing property.
The self-healing mechanism based on DA reaction could be inferred from the macroscopic and microscopic observation and molecular motion of the self-healing process. When the damaged sample was treated at 120°C, the large molecules bound by DA bonds in the sample underwent r-DA reaction and were then decomposed into  short-chain and small molecules. These short-chain and small molecules were easy to perform thermal movement and migrate so that the cracks were gradually filled by molecular chains. And then the sample was processed at 70°C, and the broken DA bonds recombined to form macromolecules, thus filling and repairing cracks.

Simulation verification
We selected three self-healing systems in table 1 to experimentally verify the MD simulation results. The conversion degrees α of the self-healing systems calculated by equation (1) from the heat flow curves (see figure 6) was in the range 86%-94%, which was similar to the conversion limit selected for the simulations, so they could be directly compared.
In the simulation, the density points corresponding to different temperatures and the curves fitted through these points are shown in figure 8. The cooling rate has a great influence on the calculation of density. The simulated cooling rate was in the order of 10 K/60 ps (1×10 13 K min −1 ), which was virtually 12 order of magnitude greater than the expected cooling rate in the experiment (10 K min −1 ). T −1 h −1 u −1 s −1 , a ratedependence study was conducted to correct the simulated Tg. The general correction is to subtract 3 K to 5 K from the difference in the predicted value every order of magnitude [35][36][37]. Therefore, the Tg of each system was corrected by reducing the MD simulations values by 36-60 K. Since the Tg value of each system does not exceed 450 K, the absolute value is small. At the same time, the cooling interval was only 500 K to 300 K, and the    density was calculated once every 25 K dropped, so the timescale of the cooling simulation was short. Under the same cooling rate, the shorter the timescale, the smaller the deviation [38,39]. As a result, we chose 36 K as the deviation value for correction. In the experiment, the variation curves of storage modulus with temperature are shown in figure 9.  The stress-strain curve can well represent the stretch properties of epoxy resin. The stress-strain curves for the mechanical and self-healing properties in the experiment and MD simulation are shown in figure 10. Comparing the results of simulations and experiments, it was found that the maximum error of four key material properties, Tg, E, UTS, and h did not exceed 4%, 14%, 15%, and 3%, respectively, as shown in figure 11. The MD simulation could replace the experiment for designing resin systems. Therefore, MD simulation results could be used as raw data set for ML.

Machine learning prediction
Through BP neural network training, the loss value of 1-100 different iterations was acquired. It can be found from figure 12 that when the number of iterations was around 60, the loss value declined to a minimum and gradually converged, which proved that at this time the error between the predicted value and the target value was the smallest and the neural network model was the most accurate. Therefore, the model after 60 iterations was selected as the well-trained neural network model to check the test data set.
As shown in figure 13, comparing the predicted values with the actual values through 5 groups of test data sets, we found that the mean absolute percentage error (MAPE) could be controlled within 8%. Thus, the welltrained neural network model is reliable and reasonable enough to predict more unknown self-healing epoxy resin systems.
The effect of composition mass ratio on performance was obtained by predicting more than 1,000 selfhealing epoxy system data sets, as shown in figure 14. In terms of Tg, E, and UTS, it can be found from the heat map that the area near the angle bisector of DGEBA and FGE presented a larger value range. The more the composition deviated from this area, the smaller the corresponding value. The value was the maximum at the bottom of the angle bisector. To keep Tg at a high level, the appropriate mass ranges per 100g mixture were 37%-100% DGEBA, 0%-28% FGE, and 0%-35% BMI, which could reach 120°C. When DGEBA was 25%-100%, FGE was 0%-37%, and BMI was 0%-38%, E of this mass ratio would increase to1.58 GPa. To obtain UTS of 27 MPa or above, the region was located on 29%-100% DGEBA, 0%-32% FGE, and 0%-39% BMI. With regard to h, it can be found from the heat map that the area near the vertical bisector of the BMI presented a larger value range. The self-healing effect in other areas is poor. h was the maximum at the top of the vertical bisector. The specific mass ratio with 0%-48% DGEBA, 24%-80% FGE, and 23%-77% BMI could obtain a large h of up to 65%.
An appropriate mass ratio between FGE and BMI is required for maximum self-healing effect. When the mass ratio of FGE and BMI was about 10/14, h was the highest. This is mainly based on the DA reaction ratio of FGE and BMI. On this basis, the mutual excess within a certain range will promote the DA reaction efficiency and increase the density of DA bond. However, if the amount was too large, the entire resin system would be damaged, which would not only degrade self-healing property but also severely damage thermal and mechanical properties. Under the condition of the same FGE/BMI, DGEBA could significantly improve Tg, E and UTS, mainly because the addition of DGEBA greatly increased the crosslinking density of the resin system, forming a dense network structure and improving its mechanical property. But at the same time, increased epoxy density reduced the density of DA bond, thus degrading the self-healing property. Therefore, an appropriate ratio can ensure self-healing property while improving mechanical property.
We comprehensively considered Tg, E, UTS and h to design a self-healing epoxy resin system in which a mass ratio of DGEBA/FGE/BMI was 39/27/34, the Tg of the designed epoxy resin was 128°C, E was 2.13 GPa, UTS was 34 MPa and h was 72%. The properties of the designed self-healing epoxy resin system were experimentally validated, and the value of each property is listed in table 3.

Conclusions
The method of designing self-healing resin systems based on MD simulation and BP neural network was used to solve the contradiction between self-healing and mechanical properties. Through ML analysis, it was found that DGEBA could significantly improve Tg, E, and UTS, suggesting that it could simultaneously ensure the mechanical and thermal properties of the self-healing epoxy resin system. However, FGE and BMI were needed to improve self-healing property, and an appropriate proportion between FGE and BMI was required to realize high h. The self-healing epoxy resin system with the optimal property was designed, where the mass ratio of DGEBA/FGE/BMI was 39/27/34. The Tg, E, UTS and h of the designed epoxy resin were 128°C, 2.13 GPa, 34 MPa, and 72% respectively.