The surface modi ﬁ cation e ﬀ ect on the interfacial properties of glass ﬁ ber - reinforced epoxy: A molecular dynamics study

: In this work, the e ﬀ ect of common functional groups, namely hydroxyl, formyl, carboxyl, and amine groups on the interfacial behavior of surface - modi ﬁ ed glass ﬁ ber - reinforced epoxy is investigated at molecular scale. The interfacial properties of the epoxy/silica coated with di ﬀ erent functional group systems are quanti ﬁ ed by performing pulling test using the steered molecular dynamics simulations. It is found that the system with hydroxyl groups has a relatively lower interfacial inter -action, exhibiting an adhesive failure mode. When partial hydroxyl groups are replaced by carboxyl, amine, and formyl groups, respectively, the interfacial interactions are increased and these systems exhibit a cohesive failure mode where failure happens in the epoxy close to inter face. A relatively higher force is required for the adhesive debonding, while more energy can be dissipated for the cohesive debonding. Because the increased interfacial interactions can prevent the mobility of polymer chains, and delay the propagation of micropores in the matrix, leading to the epoxy matrix with a high ability of energy absorption. Our work provides an insight into how func tional groups a ﬀ ect the interface debonding behavior of glass ﬁ ber - reinforced epoxy, o ﬀ ering a guideline for control of the interfacial properties of such composites through surface modi ﬁ cation techniques.


Introduction
Glass fiber-reinforced (GFR) epoxy resin has been widely used as insulating materials, adhesives, electronic packing materials, and matrix for functional composites because of its thermal stability, lightweight, superior mechanical properties, excellent resistance to corrosion, and outstanding electrical insulation [1,2]. GFR epoxy composites have revealed significant applications in a variety of engineering fields such as automobile, automotive, aerospace, and construction [3]. However, glass fiber has poor adhesion with epoxy matrix. It has been found that most failure cases of GFR epoxy composites are correlated to the weak interfacial bonding between glass fiber and epoxy matrix [4,5]. The durability of reinforced epoxy composites is highly dependent on the interface interaction. Therefore, an appropriate control over the interfacial characteristics is required in order to achieve the optimum performance of GFR epoxy composites for their safe and reliable applications.
The surface of glass fiber, modified by using alkali treatment, acetylation, electroplating, plasma treatment, and grafting, has been widely adopted to improve the interfacial adhesion between glass fiber and epoxy matrix [6,7]. For example, the surface polarity of glass fiber can be improved by plasma treatment, where large numbers of functional groups are introduced leading to an improvement in the interfacial bonding [8]. Numerous direct and indirect experimental methods have been developed to quantify the interfacial properties of fiber-reinforced composites. The most common method to investigate the interfacial properties is through fiber pull-out test, which consists of debonding, bridging, and followed pull-out process, occurring on atomic, microscopic, and macroscopic levels at the interface [9,10]. The knowledge of the sequences of events occurring on these different levels is extremely important to understand the nature of the interfacial properties. Moreover, the interface region that controls the stress transfer between fiber and matrix is primarily dependent on the level of interfacial adhesion [11]. The interfacial adhesion is confined to a several hundred nanometers-wide boundary region at the interface of fibers and polymer matrix, known as the interphase [12]. The failure of GFR epoxy initiates from the failure of interphase. The intrinsic of the surface modification of glass fiber is to modify the failure behavior of interphase. However, due to the limitation of experimental characterization, how the functional groups affect the failure of interphase is still ambiguous, making it difficult to control the interfacial properties of GFR composites. The effect of functional groups on the interfacial debonding that occurred at the early stage of fiber pull-out is still unknown, and the underlying reason for the variation in failure mechanism for surface-modified GFR composites is unrevealed.
Molecular dynamics (MD) simulations can depict the microstructure evolution, including the rearrangement of atoms, the change in polymer chain conformation, and the intra-and inter-molecular interactions [13,14]. MD simulations have become a powerful approach for exploring the dynamical processes of conformational changes and predicting mechanical properties [15]. The relationship between interfacial structures and properties can be revealed by MD simulation approach. For example, the optimal chemical functionalization groups grafting on carbon fiber to improve the mechanical properties of carbon fiber-reinforced polypropylene have been successfully determined by MD simulations [16]. The interface adsorption mechanism of surface-modified carbon fiber-reinforced epoxy has also been validated by MD simulations [17]. It has been found that oxygen atoms on the surface of the epoxy matrix are first accumulated on the carbon fiber surface; the molecular chain of epoxy is then driven towards the carbon fiber surface with functional groups [17]. The details of how functional groups affect the interfacial bonding between fiber and polymer matrix can be figured out by MD simulations.
The objective of this work is to investigate the effect of functional groups on the interfacial behavior between glass fiber and epoxy matrix at molecular scale, quantifying their effects on interfacial failure. As silica is the major composite, accounting for more than half of the weight fraction in commercial S-glass fiber, the silica structure is employed to represent glass fiber in this work. The scope of this work is to first model the initial structures of the cross-linked epoxy and modified silica with different functional groups introduced on the surface. The physical properties of the cross-linked epoxy and silica are predicted and compared with experimental results so as to validate the reliability of the modeled structures and selected forcefield. Subsequently, the bilayer structure is composed based on the equilibrated cross-linked epoxy and silica structures. The interfacial properties of epoxy/silica functionalized with different functional groups are studied by performing the pulling test using the steered MD simulations. The interfacial properties such as the force, energy, and displacement for interface failure are quantified. Finally, the conformational change of interface is captured to figure out the location of interface failure. The reason for the variation in interfacial properties of epoxy/silica systems with different functional groups is figured out; and the mechanism of interface failure for surface-modified GFR epoxy composites is discussed. Such understanding of the interactions between glass fiber with different functional groups and epoxy matrix at atomistic scale can provide valuable theoretical support for interface control of GFR epoxy composites.

Computational details
The atomistic simulations start with the construction of full-atomistic models of cross-linked epoxy and silica with different functional groups in Accelrys Materials Studio [18]. The interfacial models are constructed using Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) [19]. The equilibration and dynamic deformation are carried out using the parallel MD code LAMMPS. The details of the model construction and simulation procedures are presented.
The critical step in MD simulations is the selection of an appropriate forcefield because the forcefield determines the accuracy of the predicting properties related to atom interactions of materials [20]. The conventional polymer consistentforcefield (PCFF) is chosen. The potential energy of PCFF comprises a set of covalent-related interactions such as the bond interaction between pairs of bonded atoms, the angle interaction between three consecutive bonded atoms, the dihedral interaction, and the improper interaction between quadruplets of atoms and non-bonded interactions. Specifically, the covalent-related potential energy depends on the bond lengths, bond angles, torsion angles, and improper out-of-plane angles. The non-bonded potential energy is determined by the van der Waals interactions and Coulomb interactions. These interactions dominate the mechanical behaviors of the molecular materials. The PCFF has been successfully used to predict the structural, conformational, and vibrational properties of a broad range of molecules in condensed phases, and it also has an experimentally comparable precision in predicting molecular properties in condensed phases [21]. This potential has been successfully applied in the simulations of a wide range of organic and inorganic materials including silica and epoxy resin. In this work, PCFF is used to cross-link epoxy resin and to describe the interactions within the cross-linked epoxy and silica, as well as between epoxy and silica. Although some covalent bonds can be formed between functional groups on glass fiber and epoxy, the non-bonded interactions make a dominant contribution to the interfacial adhesion. Because the introduced functional groups on the surface glass fiber by surface modification method is mainly to improve the wettability between glass fiber and epoxy. Thus, the interfacial interaction between epoxy and silica governed by the nonbonded interactions is considered in this work.
The default high performance epoxy resin, diglycidyl ether of bisphenol A (DGEBA), is used in a wide variety of applications, including as the matrix phase in composite systems. When cured with an aromatic diamine, such as diaminodiphenyl sulfone (DDS), the properties of the resultant cured epoxy resins are much sought after, including an improved glass transition temperature and excellent chemical resistance [22]. In this study, the highperformance epoxy resin DGEBA and curing agent 4,4′diaminodiphenyl sulfone (44DDS) are selected as the representative. The β-cristobalite is served as the starting conformation for the creation of silica. This is because silica structure is constructed via the melt-quench technique and the β-cristobalite can undergo a direct phase transition to the liquid state as the temperature rises higher than the melting point [23]. The silica structure is obtained by heating the crystalline β-cristobalite to high temperature for melt and then cooling down to room temperature [24].

Cross-linking of epoxy
The molecular structures of DGEBA and 44DDS are shown in Figure 1. DGEBA reacts with the curing agent 44DDS to form a cross-linked structure. The details of the polymerization mechanism are shown in Figure 2(a). The epoxide groups in DGEBA molecules can be activated and the C-O-C bonds in the epoxide group are broken with the formation of a reactive CH 2 site [25]. The activated CH 2 site reacts with primary amine hydrogens in 44DDS, resulting in the formation of secondary amine [26]. The secondary amine in turn reacts with open epoxide groups. A crosslinked network is formed through the reaction between open epoxide groups and the amine groups.
We start to construct the cross-linked epoxy structure by packing activated DGEBA (C-O-C bonds in epoxide groups are broken and the hydroxyl groups are attached at the activated CH 2 site.) with 44DDS using the amorphous cell module in Materials Studio. The stoichiometric mixing ratio of DGEBA to 44DDS is 2:1 with 400 monomers of DGEBA and 200 monomers of 44DDS. A cubic primitive cell with a length of 6.4 nm is constructed. The energy and geometry of the system are minimized by the conjugate gradient method and then equilibrated at 300 K in the canonical (NVT) ensemble for 1 ns followed by another 1 ns equilibration in the isothermal-isobaric (NPT) ensemble at 300 K and 1 atm. Subsequently, the equilibrated system is heated to 400 K for cross-link. The condensation reaction occurs between the hydroxyl groups in open epoxide groups of DGEBA and hydrogens in amine groups of 44DDS at 400 K. The cross-linked structure is obtained with a cross-linking degree of 85%, as shown in Figure 2(b). Next the condensation reaction occurs within DGEBA monomers where the unreacted open epoxide groups change to inactivated epoxide groups through the condensation reaction between hydroxyl groups, meaning that DGEBA molecules are ended with epoxide groups. The system is finally cooled down to 300 K for equilibration in NPT ensemble. After equilibration, the obtained cross-linked structure is served as the initial structure of epoxy resin with a size of 6.4 nm × 6.4 nm × 6.4 nm.
The physical properties of epoxy resin are predicted in LAMMPS. The cross-linked structure obtained from Materials Studio is first minimized and equilibrated for 1 ns in NVT ensemble at 300 K, then equilibrated for another 1 ns in NPT ensemble at 300 K and 1 atm. The purpose of doing NVT before NPT is based on algorithmic stability. Velocity generation in MD simulations can crash if coupled directly with a barostat. The equilibration is often done under NVT for a period to get the velocity distribution reasonable, followed by NPT. The root-meansquare displacement (RMSD) is checked; it is found that the value of RMSD at the end of 500 ps is almost constant indicating that the system is equilibrated. The density of the cross-linked epoxy is about 1.14 g/cm 3 , located in the range of 1.1-1.4 g/cm 3 tested by experiments [25,27]. The glass transition temperature (T g ) and Young's modulus (E) of the cross-linked epoxy are predicted. Specifically, the equilibrated structure is heated to 690 K and then cooled down to 200 K at a temperature step of 30 K with a cooling rate of 0.5 K/ps. At each temperature, an equilibration of 500 ps is performed using NPT ensemble at 1 atm, so that the predicted results can be directly compared with experimental data. The specific volume as a function of temperature is obtained to evaluate T g of cross-linked epoxy. The tensile deformation of the system is performed at 300 K and 1 atm where the system is deformed along z direction with a fixed strain rate of 10 8 /s. The system is equilibrated for about 50 ps following each step of deformation. The stress under a specific strain is obtained from the average stress for the last 5 ps equilibration. The value of E is predicted from the stress-strain curve.

Construction of silica model
The β-cristobalite has a lattice constant of 7.16 Å with an oxygen-to-silicon ratio of 2:1 shown in Figure 3(a). The crystalline β-cristobalite system is cleaved on (001) plane. The dangling oxygen atoms on the (001) surface are all saturated by the hydrogen atoms and dangling silicon atoms are saturated by hydroxyl groups, as shown in Figure 3(b). The density of the surface hydroxyl groups is about 7.9/nm 2 , located in the range of the experimental results [28]. It is found that the surface of silica is covered, to some extent, with hydroxyl groups during the formation processes; and the overall coverage of hydroxyl groups on a surface ranges from 0 to 9.4/nm 2 in experimental tests [29]. The silica is first equilibrated at 300 K in NVT ensemble for 1 ns, then equilibrated in NPT ensemble at 300 K and 1 atm for another 1 ns with the periodic boundary condition. The equilibrated system is then heated to 5,000 K for melting and equilibrated for 1 ns in NPT ensemble. Next the equilibrated silica is cooled down to 300 K at a rate of 5 K/ps and equilibrated at 300 K at 1 atm for 1 ns. The density of constructed silica is about 2.18 g/cm 3 , close to the experimental data [30].
In order to study the effect of functional groups on the interfacial bonding, half of the number of hydroxyl groups on the surface of silica are randomly substituted by functional groups such as formyl groups, carboxyl groups, and amine groups with a density of 4/nm 2 . The silica models with different functional groups are equilibrated at 300 K and 1 atm in NPT ensemble for 1 ns. The RMSD is checked, and it is constant at the end of 500 ps, which indicates that the system is in an equilibrated state.
Radial distribution function (RDF) analysis is among one of the most important methods for revealing the structural features of a system and is calculated to quantify bond types. RDF is the probability density of finding atoms A and B at a distance r averaged over the equilibrium trajectory. Figure 4 shows the partial RDF of different silica systems. In the silica system with hydroxyl groups, the Si-O bond in inner silicon oxygen tetrahedron (bulk Si-O) and in silanol groups (Si-OH) are about 1.58 Å, as shown in Figure 4(a). This has a good agreement with the value (about 1.61 Å) measured by X-ray diffraction [31,32]. The length of bulk Si-O bond in silica models covered with formyl, carboxyl, and amine groups are similar to that in the silica models covered with hydroxyl groups as shown in Figure 4(b). The bond length of Si-N in the silica model covered with amine groups is around 1.68 Å, consistent with the value of 1.69 Å measured by experiments [33]. The length of Si-C in the silica model covered with formyl and carboxyl groups are 1.98 and 1.93 Å, respectively, in Figure 4(c), similar to those tested by experimental and computational approaches [34].

Construction of epoxy/silica structure
The equilibrated epoxy, silica with different functional groups on the surface, and a vacuum layer with a thickness of 2 nm are sequentially placed in order from bottom to top with the formation of the complete simulation systems. The interlayer separation between silica and epoxy is initially selected to be 5 Å, which is subsequently adjusted by the repulsive and attractive forces at the interface during equilibration process. The periodic boundary conditions are applied in x, y, and z directions to avoid the finite-size effects. The vacuum layer added is in order to avoid the interfacial interaction caused by the periodic boundary condition in z direction. The interactions between epoxy and silica mainly consist of non-bonded interactions including the Coulomb interaction and van der Waals interaction. The constructed epoxy/silica system is first energy minimized, and then equilibrated at 300 K in NVT ensemble for 1 ns, followed by another 1 ns equilibration at 300 K and 1 atom in NPT ensemble. Finally, the RMSD is checked, and it is constant at the end of 500 ps, indicating that the system is in an equilibrated state.
The interfacial bonding energy of bilayer materials is calculated based on the formula [26]: where E bonding is the interfacial bonding energy of the epoxy/silica system; E epoxy is the potential energy of the equilibrated epoxy system; E silica is the potential energy of the equilibrated silica system; E bilayer is the potential energy of the equilibrated epoxy/silica system and A is the interface area. The interfacial debonding behavior of different bilayer systems is studied by pulling test using steered molecular dynamics (SMD) simulations. The SMD simulations, based on the principles of atomic force microscopy technique, can provide the conformational change in atomistic scale and energy variation during dynamic deformation [35]. In SMD simulations, the center-of-mass (COM) for epoxy and silica is attached by a virtual spring. The atoms are displaced by applying a constant velocity (v) to them. A restoring force is applied to the atoms and the magnitude of the force is related to the spring constant (k). The epoxy and silica are debonded along z direction under equivalent opposite force perpendicular to the interface. The virtual spring force is determined by ref. [36]: where t is the time; n is a unit vector for the direction of pulling; R t and R 0 are the displacement at t and the initial displacement between the COM of epoxy and silica, respectively. U and F are the potential energy and virtual spring force, respectively. As shown in equations (2) and (3), the pulling force applied is dependent on the spring constant and the pulling velocity. If the spring constant is too small, the interface cannot be debonded. Otherwise, if the spring constant is too large or the pulling velocity is too fast, there is a problem to sample the interaction between fiber and epoxy along the pulling displacement.
On the other hand, the velocity cannot be too small, considering the simulation time. Taking all these into consideration, we use the velocity of 1 Å/ps and a spring constant of 100 kcal/(mol · Å 2 ) to ensure the high simulation accuracy and reasonable computational cost of simulations. The epoxy and silica are debonded along z direction under equivalent opposite force perpendicular to the interface. The pulling force F is collected every 0.1 fs. The initial length of the spring is equal to the initial displacement between the COM of epoxy and silica. As the value of COM for epoxy and silica along x and y directions are nearly the same, the variation in the distance between the COM of epoxy and silica along the z direction can be regarded as pulling displacement (d). This means that there is relationship of ( ) . The free energy change is equal to the work done on the system in a reversible isothermal process and equal to the integral of an externally applied force over the coordinate. During SMD simulations, the free energy change is defined as the potential of mean force (PMF). The SMD simulations are performed in NVT ensemble at 300 K. Because SMD simulations require an empty volume along one Cartesian axis, the system has free space to stretch without colliding with its periodic images [37]. To avoid this problem, a vacuum slab is added, and the NVT ensemble is used for SMD simulations to avoid a barostat control that can reduce the empty space added. The conformational changes in interface are captured with OVITO. Figure 5 shows the partial RDF of bonds in epoxide groups and the cross-linked bonds between DGEBA and 44DDS. During the cross-linking process, the C-O bond in the epoxide group is broken as shown in Figure 5(a).

Physical properties of cross-linked epoxy
The length of C-C and C-O bonds is 1.48 and 1.43 Å, respectively, as shown in Figure 5(b). The C-N bond is formed between DGEBA and 44DDS for cross-linking; its length is about 1.53 Å as shown in Figure 5 The physical and mechanical properties of the crosslinked epoxy have been predicted and compared with the experimental data for validation. Figure 6(a) shows the temperature dependence of the specific volume for the cross-linked epoxy during the cooling process. The specific volume of the cross-linked epoxy at each temperature has been calculated from the average density of the system for the last 50 ps of equilibration. The abrupt change in the slope of the curve determines T g . The intersect of the extrapolations denotes that T g of the cross-linked epoxy is about 489 K, located in the range from 460 to 530 K tested by the experimental approach [22,26,39]. Figure 6(b) shows the stress-strain curve of the tensile deformation for the cross-linked epoxy. Young's modulus is about 4.4 GPa. This is located in the range of experimental results and close to other computational results [40,41]. All these indicate the reliability of PCFF in terms of predicting the mechanical properties of the modeled cross-linked epoxy. As the functional groups introduced on the surface of glass fiber mainly affect the conformational change in epoxy during debonding, a reasonable selection of forcefield to describe the epoxy properties is important.

Interfacial properties of epoxy/silica systems
The effect of interfacial properties between epoxy and functionalized silica is investigated by performing the pulling test. The pulling force and PMF as a function of pulling distance are shown in Figure 7. In the pulling test, the force first increases linearly with pulling displacement, and the functional groups effect can be negligible as shown in Figure 7(b). The force, then, increases   . The pulling displacement for the different systems at the point of F max is similar, around 1.1 Å. After pulling force reaches the maximum, it reduces significantly with the increment in pulling displacement. The value of PMF as a function of d is shown in Figure 7(c). PMF first increases linearly, and then increases progressively to reach the first peak; PMF finally keeps increasing with fluctuation. It is found that when PMF reaches the first peak, the pulling force is almost reduced to zero. . The change in interfacial interaction (E inter ) between epoxy and silica with different functional groups during debonding is shown in Figure 7(d). It is clear that the interfacial interaction is increased when partial hydroxyl groups on the surface of silica are replaced by formyl, carboxyl, and amine groups. There is a relationship of . The E inter reduces and then reaches a stable value from the point of PMF peak . This indicates that the bilayer system is completely debonded at the point of PMF peak . Moreover, when the interface is completely debonded, the interaction between epoxy and silica with hydroxyl groups is reduced to zero, while that between epoxy and silica with formyl, carboxyl, and amine groups still exists. This indicates that the surface modification of glass fiber changes the interface failure mode.
The interfacial bonding energy, the value of F max , PMF peak , and corresponding d 0 for different bilayer systems are shown in Table 1 . The relationship of the interfacial bonding energy in different systems is consistent with the relationship of F max . The epoxy/silica with hydroxyl groups system has the highest force for debonding, but the lowest absorbed energy. For the surface-modified silica, the variation trend of maximum pulling force is approximately consistent with that of interfacial bonding energy and PMF. For example, compared with the system with carboxyl groups, the system with amine groups requires higher force and energy for interface debonding. In order to understand the reason for the difference, the evolution of interfacial structures during pulling is captured.
The conformational changes at different displacements are represented to figure out the reason for the change in force and energy. Figure 8 shows the conformation of the epoxy/silica system with hydroxyl groups system at different pulling displacements. After the force reaches the peak, with the continuous increment in displacement, micropores are generated between epoxy and silica as shown in Figure 8(b). These micropores grow (Figure 8c) resulting in the debonding between epoxy and silica at the point of PMF peak (Figure 8d). The interface fails at the contact surface between epoxy and silica; the debonded interface is relatively smooth as shown in Figure 8(e). The conformational change in the epoxy/ silica with formyl, carboxyl, and amine groups during pulling is shown in Figure 9. Different from the system with hydroxyl groups, the interface fails in the epoxy close to the interface for these surface-modified systems. Specifically, the interface failure of the system with formyl groups is mainly due to the localized growth of micropores as shown in Figure 9(b). When the interface fails, some polymer chains of epoxy are still bonded to the silica as shown in Figure 9(c). For the system with carboxyl groups, the growth of micropores is associated with more micropores formed (Figure 9e). The number of micropores is much larger than that in the system with formyl groups at the same displacement. The coupled effect of formation and growth of micropores leads to the failure of interface with a few polymer chains bonded to silica (Figure 9f). However, after the micropores formed in the system with amine groups (Figure 9g), the size of some micropores is reduced with the increment in displacement as shown in Figure 9(h). A great number of polymer chains are still bonded to silica when the interface fails as shown in Figure 9(i). This indicates that a stronger interfacial interaction can prevent the mobility of polymer chains, inhibit the generation and growth of micropores, and even repair the micropores.
For a better understanding of the conformational change in epoxy during deformation, the radius of gyration (R g ) of the epoxy is quantified and calculated by: The interfacial properties of different epoxy/silica systems, including the functional groups density, interfacial bonding energy (E bonding ), the maximum pulling force (F max ), the PMF for debonding (PMF peak ), and the corresponding pulling displacement (d 0 )  where M is the total mass of epoxy; r cm is the COM position of the group; r i is the position of atom i; m i is the mass of atom i; and the sum is the overall atoms in the epoxy. This parameter can be used to evaluate the degree of compression and stretch for polymer chains. Figure 10 shows the curves of R g as a function of pulling displacement d. The R g of the system with hydroxyl groups first increases, and then almost keeps constant as shown in Figure 10(a). However, it reduces significantly to reach the minimum value. This rapid reduction is correlated to the failure mode. When the interface between epoxy and silica is debonded, the motion of polymer chains in epoxy is no longer constrained by the interfacial interactions, and R g reaches the minimum at the point of PMF peak . The R g of the system with formyl groups is first increased followed by a decrease as shown in Figure 10(b). The reduction of R g is correlated to the constraint of polymer chains from the silica is reduced with the growth of micropores. The micropores generated in the epoxy close to the interface cause the increment in R g . When the interface fails, R g reaches to the minima due to a significant release of constraint. For the system with carboxyl groups, the R g fluctuates obviously before the interface fails in Figure 10(c). This is correlated to the formation and growth of multi micropores at the contact surface between epoxy and silica and in the epoxy near the interface. The R g of the system with amine groups first increases and then reduces slightly in Figure 10(d). The reduction is correlated to the formation of micropores. However, different from the system with hydroxyl, formyl, and carboxyl groups, the R g keeps growing with the increase in the displacement after the point of PMF peak . This is mainly correlated to the increased number of polymer chains bonded to silica, where there is still a strong interfacial interaction between epoxy and silica when the interface is completely debonded. Generally, R g first increases rapidly when the force grows linearly at the beginning stage. The reduction in R g is correlated to the constraint release from silica caused by the generation and growth of micropores near the interface. However, a strong interfacial interaction promotes the growth of micropores within the epoxy matrix far from the interface, leading to an increment in R g .

Failure mechanism of surface-modified GFR epoxy
The functional groups introduced on the surface of glass fiber change the interfacial interaction, significantly affecting the interface debonding behavior. As shown in Figure 11, there are two types of debonding failure mode for GFR epoxy at nanoscale: adhesive failure at the interface between epoxy and glass fiber, and cohesive failure in the epoxy matrix. Generally, higher stress is required for the adhesive failure mode, while more energy is absorbed for the cohesive failure mode. Because the increased interfacial interactions between epoxy and silica can prevent the debonding that occurred at the contact surface between epoxy and silica, the matrix has the lowest strength, causing a cohesive failure mode. A stronger interfacial adhesion prevents the mobility of polymer chains, inhibiting the generation and growth of micropores near the interface, leading to more micropores generated in the matrix far from the interface. More energy can be absorbed for the interface debonding. For example, the introduction of amine groups on the surface contributes to the improvement in the ability of energy absorption, resulting in a much larger displacement for interface debonding. Such failure mode is expected because GFR epoxy composites have a high ability of energy absorption and the components made up of such composites can avoid the catastrophic failures caused by the localized Figure 10: R g of epoxy as a function of displacement during pulling for (a) the system with hydroxyl groups; (b) the system with formyl groups; (c) the system with carboxyl groups; and (d) the system with amine groups. R g increases obviously at the beginning stage of pulling, consistent with the linear increment in pulling force. When the interface fails, R g reaches the local minima.
failure, making them reliable. Moreover, the cohesive failure of interface debonding can make it difficult for fiber pull-out that initiates from debonding between epoxy and glass fiber.
Additionally, when the interfacial debonding occurs in the epoxy close to the interface, methods can be applied to improve the properties of epoxy. For example, the crosslinking degree of the epoxy matrix can be improved to decrease the speed of micropore growth. The reinforcements such as carbon nanotubes [42] and graphene [43] can be added to inhibit micropore growth in the matrix for optimizing the performance of GFR composites.
In this work, we mainly consider the non-bonded interactions between glass fiber and epoxy matrix. In practice, the glass fiber surface modification changes not only the interfacial interaction but also the surface morphology of fiber and chemical reactivity along the fiber-matrix interface [44]. The surface roughness and the chemical bonds between epoxy and surface-modified glass fiber are not considered, because the main purpose of modification of the surface characteristics of glass fiber is to improve the compatibility or wettability between glass fiber and epoxy matrix. The naturally hydrophilic epoxy is not inherently compatible with hydrophobic glass fiber. Composites with improved wettability can have superior interfacial bonding. Moreover, the effect of roughness on interfacial adhesion is complicated. An increment in surface roughness of glass fiber contributes to the mechanical interlocking between fiber and polymer matrix. However, with a further increase in surface roughness, the fiber is damaged with a reduction in mechanical properties. Compared with the non-bonded interfacial interactions, the chemical bonds at interface and the surface roughness are more difficult to be controlled. For surface modification techniques, the factors affecting the interfacial bonding between fiber and matrix include mechanical interlocking, physical binding, and chemical bonding. Our simulated work on the physical binding effect between glass fiber and polymer matrix lay a foundation to realize precise control over the interface adhesion of GFR composites. For the example of surface modification technique of plasma treatment, our results provide a guideline for the selection of the atmosphere parameter. A nitrogen atmosphere can be selected to introduce the amine groups on the silica surface for surface modification of glass fibers.

Conclusion
The surface modification techniques with the introduction of different functional groups have been widely applied to improve the interfacial bonding between the glass fiber and epoxy resin. However, as the underlying mechanism is still not fully understood, the surface modification is still based on the trial-and-error approach to obtain the optimum performance of GFR epoxy. The enhancement of interfacial properties for GFR epoxy by surface modification is limited. In this work, the effect of functional groups on the interfacial behavior between epoxy and silica at atomistic scale has been studied. The DGEBA and 44DDS are cross-linked to form epoxy resin. The interfacial properties of the epoxy/ silica systems have been predicted by the pulling test using SMD simulations. Based on the findings, several conclusions can be drawn: 1) The epoxy/silica system with hydroxyl groups requires the highest force, but the shortest pulling displacement and the lowest energy for interface debonding. The Figure 11: The surface-modified GFR epoxy composites exhibit two types of interface debonding: adhesive failure mode and cohesive failure mode. The introduction of functional groups significantly affects the interface failure mode. The strong interfacial interactions between epoxy and silica prevent the debonding at the contact surface between epoxy and silica, leading to a cohesive failure mode.
system fails at the contact surface between epoxy and silica, exhibiting an adhesive failure mode. 2) When partial hydroxyl groups are replaced by formyl, carboxyl, and amine groups, respectively, the system fails in the epoxy close to the interface, representing a cohesive failure mode. Although the pulling force for debonding is reduced compared with that in the system with hydroxyl groups, more energy is absorbed for the rearrangement of polymer chains in epoxy. 3) When the epoxy/surface-modified silica system exhibits a cohesive failure mode, the force and energy for debonding are increased with the increment in the interface interactions. Compared with the system with carboxyl groups, the systems with amine and formyl groups require a relatively higher force and energy for debonding because of the stronger interfacial interaction.
Our work reveals the mechanism of different functional groups on the interfacial behavior of GFR epoxy. Such understanding is fundamental to figure out the reason for the variation in interfacial performance caused by different glass fiber surface modification techniques. The findings in this work lay a foundation for the accurate control of interfacial properties of GFR composites by surface modification.