Effect of the Nanoparticle Functionalization on the Cavitation and Crazing Process in the Polymer Nanocomposites

The nanoparticle (NP) functionalization is an effective method for enhancing their compatibility with polymer which can influence the fracture property of the polymer nanocomposites (PNCs). This work aims to further understand the cavitation and crazing process, hoping to uncover the fracture mechanism on the molecular level. By adopting a coarse-grained molecular dynamics simulation, the fracture energy of PNCs first increases and then decreases with increasing the NP functionalization degree α while it shows a continuous increase with increasing the interaction εpA between polymer and modified beads. The bond orientation degree is first characterized which is referred to as the elongation. Meanwhile, the stress by polymer chains is gradually reduced with increasing the α or the εpA while that by NPs is enhanced. Furthermore, the percentage of stress by polymer chains first increases and then decreases with increasing the strain while that by NPs shows a contrast trend. Moreover, the number of voids is quantified which first increases and then decreases with increasing the strain which reflects their nucleation and coalescence process. The voids prefer to generate from the polymer-NP interface to the polymer matrix with increasing α or εpA. As a result, the number of voids first increases and then decreases with increasing α while it continuously declines with the εpA. In summary, our work provides a clear understanding on how the NP functionalization influences the cavitation and crazing process during the fracture process.


INTRODUCTION
It is important to reinforce the polymer nanocomposites (PNCs) for developing highly functional materials, such as fuel-efficient tires or oil-sealing rubber parts. [1,2] Generally, addition of a large amount of nanoparticles (NPs) into the matrix can reinforce the PNCs which however increases the difficulty in uncovering the physical mechanisms. Many parameters influence the reinforcement of PNCs such as the shape, size, content and dispersion of NPs, the polymer-NP interaction and so on. [3−7] The evolution process of the cavitation and crazing plays a key role in determining the mechanical properties which helps to uncover the fracture mechanism. However, a fundamental understanding on a molecular scale has not been achieved yet.
Crazing is a unique mode of failure for PNCs, during which a strong dilatational component induces a rapid volume gain before failure. Nanovoids generate, grow and coalesce into a microcavity during the fracture process. Finally, the percolated voids induce the failure of PNCs. The slippage and extension process of polymer chains and NPs requires significant energy, making it a source of the toughness of PNCs that is important to their use in load-bearing applications. Thus, the nucleation and growth of nanovoids are vital to understanding the fracture mechanism. Recently, lots of works have focused on this fundamental question with the aim of realizing better performance. By examining the propagating crazes, the developed model predicts that capillary forces drive the shape and distribution of forming fibrils, from which one would expect surface tension to be the dominant contributor to stress. [8] The fracture toughness of glassy polymers has been investigated via constitutive equations to understand the competition between shear yielding and crazing. The development of crazes is found to depend upon the kinetics of local plastic deformation. [9] By analyzing the craze formation and propagation in glassy PNCs, three stages of the rearrangement process of fillers are revealed. [10] In addition, the remalleable, healable, and biodegradable polymeric materials with integrated high strength and high toughness are obtained which is attributed to the interpenetrating threedimensional supramolecular clusters. [11] By adopting the small-angle X-ray scattering (SAXS), the formation of voids was detected in the carbon black filled styrene-butadiene rubber. [12−15] However, the only voids with size from 20 nm to 30 nm can be detected, which is larger than that (<10 nm) of the initial ones. Furthermore, the experimental instrument is difficult to record the nucleation and the evolution process of voids. Molecular dynamics (MD) simulation offers a unique approach to characterize the microstructural evolution during the fracture process. By adopting MD simulation, the nucleation of voids prefers to appear in the region of high chain ends density [16] or the low elastic modulus. [17] In addition, the amounts of the dissipated energy can be roughly divided into three types of rapid motion: cavitation, plastic yield and bridge rupture. [18] By investigating the early behavior of pure polymer crazes, the characteristic length of stretched chains is one-third the entanglement length and the drawing stress of the system is carried primarily by covalent bonds. [19] Meanwhile, the Voronoi volume can predict the void formation and the failure in rods filled PNCs. [20] Then, nanorods incorporated into the crazes rapidly orient themselves to match the direction of the polymer fibrils while those in bulk regions remain randomly oriented. Furthermore, there is a strong correlation between local elastic moduli and the type of nonaffine plastic failure that is otherwise associated with polymer glass deformation. [21] NPs can act as temporary cross-link bonds which thus form the polymer-NP networks to enhance their mechanical toughness. [22] It is reported that small fillers with the high mobility can dissipate the deformation energy via an improved release of local tension. [23] However, the reinforcement is nearly independent of the NP mass which will affect their mobility. [24] The optimal fracture property is realized at the moderate grafting density which can be explained by analyzing the stress contribution of different components. [25] NPs tend to form the contact aggregation structure which is the main obstacle for reinforcement. Among the various methods, [26−28] the chemical functionalization of NPs seems much more effective than others to improve their dispersion and reinforce the polymer. However, a fundamental understanding of the cavitation and crazing process has not been clearly identified yet. To understand the effect of the NP functionalization on the fracture behavior of PNCs, it is vital to characterize the nucleation and evolution process of voids. In this work, a coarse-grained model is used to investigate their fracture behavior at the molecular level. The structural evolution is characterized in details with respect to the strain, which establishes their relationship with the stress-strain behavior. The following two questions are hopefully to be answered: (1) How the NP functionalization degree and the interaction between polymer and the modified beads affect the stress contribution by polymer chains and NPs to the total stress respectively. (2) Where the nucleation of voids occurs and what the evolution process of voids is.

MODELS AND SIMULATION DETAILS
The coarse-grained model of the functionalized spherical NPs filled PNCs is selected in this work. The classic bead-spring α ε σ model [29] is used to simulate polymer chains. Each chain contains thirty beads whose length is less than the entanglement length. There are 1000 chains for each system. The NP is modeled as the spherical shell with 48 Lennard-Jones interaction beads distributed evenly on the surface. There are two types of randomly distributed beads on the functionalized NP, which are the strong (A) and the weak (B) type. The functionalization degree of NPs ( ) is defined as the ratio of the number of A beads to that of total beads in one NP (48 in this work). Each system of size ≈ 34σ contains 112 NPs. It is roughly corresponding to the volume fraction of NPs = 9.22% which is defined as the ratio of the volume of NPs to the box volume. The diameter and mass of the polymer bead and the surface bead of NPs are same, which are denoted by σ and m, respectively. The diameter and mass of the inner center of NPs is 0σ and 16m, respectively, which ensures the same average density of our NPs and the homogeneous NPs with the same diameter. [30] Fig. S1 (in the electronic supplementary information, ESI) presents the snapshot of NPs for different α. When mapping the coarsegrained model to real polymers, the interaction parameter is set to be about 2.5−4.0 kJ·mol −1 for different polymers. [29,31] Meanwhile, one polymer bead with a diameter 1 is about 1 nm which roughly corresponds to 3−5 repeating units of polybutadiene. [32−35] In total, our simulation falls within a parameter range in experiments which can capture these typical polymer systems.
The truncated and shifted Lenard-Jones (LJ) potential is used to describe the non-bonded interactions between all beads [29] as follows: where r is the distance between two interaction sites. C is a constant to guarantee that the potential energy is continuous. r cutoff stands for the distance at which the interaction is truncated and shifted so that the energy is zero. The polymerpolymer interaction parameter and its cutoff distance are and . The interactions between polymer and the surface A beads are set to be strong interaction ( and ) so that A beads are miscible with polymer chains. The interaction parameter and the cutoff distance between polymer and the surface B beads are and . The surface bead-surface bead interaction parameter and its cutoff distance are and . The pair interaction energy scale and the length scale of our model are denoted by and , respectively. It is noted that our aim is not to consider a specific polymer chain. Thus, the reduced LJ units σ and ε are set to be unity, which means that all the simulated quantities are dimensionless.
The bonded interaction between the adjacent beads (including both chains and NPs) is described by the stiff harmonic potential, which is given by The are set to be 1 for polymer bonds and the surface bead-surface bead bond while it is 2.0 for the surface bead-NP core bond. These parameters can avoid high-frequency modes and chain crossing while keeping a certain stiffness of the bonds. The harmonic potential is proved to be efficient in modeling polymer systems. [36−38] After determining the force field parameters, the simulations are started from a nonoverlapped configuration of all the polymers and the NPs in a large box. [39,40] Subsequently, the NPT ensemble is used to compress the system for 5000τ where τ is the reduced time unit. The Nose-Hoover temperature thermostat and pressure barostat are adopted to fix the temperature and pressure to be and , respectively. The equations of motion are integrated by the velocity-Verlet algorithm with the time step of . Periodic boundary conditions are employed in three directions. Further equilibration is performed under the NVT ensemble with for 20000 . After that, the is gradually reduced to be 1.0 by 0.02 every 250τ. Next, the systems are equilibrated under NVT with for 20000τ. It has been checked that each chain has moved at least 2 . The number density of polymer beads is about 0.85 which corresponds to the density of melts. After that, the structure and dynamics data are collected for the ensemble average. To simulate the fracture behavior, the tri-axial deformation is adopted to induce the failure of PNCs. The length of the simulation box is extended in the tensile direction while they are fixed in other two dimensions which induces a positive effective stress in all directions. [20,41] Periodic boundary conditions are still applied in three directions during the fracture process. The strain rate is which is adopted by Gao et al. [42] It is noted that our adopted strain rate is much larger than the experimental value. Each system is deformed independently along x, y and z directions, respectively, which can give the average stressstrain curve. All simulations have been performed by using a large scale atomic/molecular massively parallel simulator (LAMMPS). [43] In addition, the nucleation and evolution process of voids is a key factor to understand the fracture process of PNCs. Thus, the method is introduced to define the nanovoids. First, the simulation box is divided into small cubic sub-cells of size Δ. Then, whether the polymer beads or nanofiller beads are within these small sub-cells will be checked, which depends on their positions. The small subcells are unoccupied which act as voids if there are no beads within them. Otherwise, they are not belonging to voids. At this moment the position and the number of sub-cells are obtained. Then, the unoccupied sub-cells are belonging to the same void if they are neighbors. The number of voids is obtained from this analysis. The size Δ of the sub-cell is about 1.7 in this work, which is comparable with that in the previously published work. [16] Meanwhile, there are no voids for the unstrained system.

RESULTS AND DISCUSSION
Nanoparticle Functionalization Degree α α The functionalization of NPs is an effective method for enhancing their compatibility with polymer chains. This will affect the dispersion state of NPs which is correlated with the fracture property. Thus, we mainly investigated the effect of the nanoparticle functionalization degree on the fracture property of PNCs in this section. Before discussing the stress-strain behavior, it is very important to characterize the dispersion state of NPs for different NP functionalization degree by calculating the radial distribution function (RDF) in Fig. S2 (in ESI). The peak reflects the direct contact aggregation of NPs while that at stands for the sandwiched structure of NPs via one polymer layer. The direct contact aggregation of NPs is formed at which can be clearly proved by the high peak at . This is induced by the depletion effect which generates from the polymer chains. Then, this depletion effect is gradually reduced with increasing the , which breaks the direct contact aggregation down. Accompanied by it, the peak at disappears at = 0.2. Interestingly, the peak at increases with further increasing α, which indicates that NPs tend to form the aggregation sandwiched via one polymer layer. The snapshots of NPs are presented for different in Fig. 1 to intuitively observe their dispersion state. The most uniform dispersion of NPs appears at the mediate . In total, NPs form the contact aggregation at , the dispersion structure at , and the bridging or tele-bridging structures at , which is consistent with the results by the microscopic polymer reference interaction site model. [44] Then, we performed the tri-axial deformation to analyze the microstructural evolution during the fracture process, which establishes their relationship with the property. First, the dependence of the stress on the strain is analyzed for different , which is shown in Fig. 2(a). It is observed that the stress first shows a linear increase with the strain and then is gradually reduced to be zero. The strain at the maximum stress is roughly 0.1. To quantitatively characterize the fracture property, we calculated the maximum stress, the elongation (the strain at stress = 0.0) and the fracture energy (area formed by stress-stain curves and the x axis), which are summarized in Table SI (in ESI). It is observed that the maximum stress gradually increases with the because of the gradual dispersion of NPs at the strong polymer-NP interaction. However, the elongation shows first an increase and then a decrease. This is because of more stress concentration on the polymer-filler interface, which induces the heterogeneous stress distribution. It restricts the chain deformation, which reduces the elongation at . As a result, the fracture energy first rises and then declines, which indicates the optimal property at the moderate = 0.4. In addition, compared with the neat system (only containing α chains), the reinforcement effect is not observed at , which is proved by the low fracture energy while it is realized at . Then, the bond orientation degree of polymer chains is characterized with the strain to further understand the fracture behavior. It is denoted by the second-order Legendre polynomials , which is defined as . Here, the denotes the angle between the bond vector and the tensile direction. The change of the with the strain is presented in Fig. 2(b) for different . The polymer bonds are randomly orientated at strain = 0.0, which is proved by the = 0.0. Then, first rises and then declines to 0.0 with increasing the strain for different . This indicates that the chains are first stretched and then contracted during the fracture process. Meanwhile, the strain at the = 0.0 is maximum at the moderate = 0.4, which is consistent with the elongation. However, the maximum is less than 0.05, which exerts a limited effect on the stress-strain behavior. Then, the number of the nearest neighbor NPs surrounding one NP at a separation closer than is calculated, which is denoted by N num . This reflects the network structure of NPs whose variation with the strain is shown in Fig. S3 (in ESI). The results present that the N num first decreases and then increases with increasing the strain for different . This indicates that the distance between NPs is first enlarged induced by the deformation. Then, accompanied by the contraction of polymer chains, their α α α α α α α distance is reduced. We next analyzed the stress contribution of polymer chains and NPs to the total stress respectively to understand how the affects the fracture property. From Fig. 3, it is observed that the stress by NPs is gradually enhanced with the increase of , which is attributed to their improved compatibility with polymer. However, the stress by polymer chains does not exhibit a continuous increase. This is because the stress by polymer chains is strongly negative at strain = 0.0 especially for = 1.0. The attractive interaction between polymer and NPs produces an organized region of polymer around the NPs. The NPs pull polymer chains inward and the shell of polymer around NPs is tightly packed. This leads to the negative stress by polymer chains and the positive stress by NPs at strain = 0.0. [45] The difference of stress by polymer chains at strain = 0.0 and 0.1 gradually rises from 0.33, 0.81, 1.01, 1.04, 1,07 to 1.11 with the . Next, the percentage of the stress by polymer chains and NPs is calculated respectively for different in Fig. 4. It is found that the nearly 100% stress is borne by polymer chains because of their weak interface at = 0.0. Then, the percentage of stress by polymer chains first increases and then decreases with the strain while that by NPs shows a contrast trend. Meanwhile, the polymer chains bear a less percentage of stress with increasing the while NPs bear a more one. In total, the contact aggregation of NPs and the weak interface lead to the low fracture energy at Last, the voids are quantified during the fracture process, which is very important to understand the fracture behavior. The change of the number of voids with the strain is presented in Fig. 5(a) for different . Meanwhile, Fig. 5(b) shows the snapshots of systems corresponding to some typical strains to intuitively observe the voids. No voids are observed at strain = 0.0. Then, some single voids are generated in the matrix with increasing the strain (< 0.10). Then, these single voids gradually grow into large voids and the number of the voids increases to the maximum value with further increasing the strain. Next, the small voids are met with each other to coalescent into a large one. The number of voids begins to decrease when the coalescent rate is larger than the generation rate of new voids. At last the percolated voids appear, which induces the failure of PNCs. In addition, the number of voids at strain ≈ 0.22 is recorded for different in Fig. 6(a) which first increases and then decreases with the . To under-≈ 1.0σ α α α α stand it, we probed the positions where the voids generate. The initial single voids are considered to generate at the polymer-NP interface when they meet the NPs within the surfacesurface distance . Otherwise, they generate in the matrix. The results are presented in Fig. 6(b). There is a high probability of voids to generate at the polymer-NP interface at = 0.0 which is a weak region. Then, the polymer-NP interface is reinforced with the increase of , which improves the probability of voids in the matrix. The region of the polymer-NP interface is limited at = 0.0, which leads to the small number of voids. The voids can be generated at the interface or in the matrix at 0.2, which results in the maximum number of voids. Last, the voids can be generated in the matrix, but not at the interface at = 1.0, which reduces the number of voids again. These can rationalize the observed results.

Interaction between Polymer and A Beads
After discussing the effect of the NP functionalization degree on the fracture properties in detail, we turned to understand the effect of the interaction between polymer and A beads in this section. The is fixed to be 0.4 where the PNCs exhibit the maximum fracture energy. The varies from 1.5 to 5.0, which simulates the different interaction strengths. Similarly, the dispersion state of NPs is analyzed for different by calculating their RDF in Fig. S4(a) (in ESI). Combined with the snapshots of NPs in Fig. S4(b) (in ESI), the peak at gradually increases with the , which indicates more bridging structure of NPs via one polymer layer. Then, the stress-strain curves are shown in Fig. 7 for different . As summarized in Table S1 (in ESI), both the maximum stress and the elongation rise with the increase of , which leads to the improved fracture energy. Similarly, the bond orientation degree of polymer chains is first analyzed as a function of the strain in Fig.  S5 (in ESI). It shows that the strain at = 0.0 increases with the , which is consistent with the elongation. Then, the stress contribution of polymer chains and NPs to the total stress is calculated respectively, which is shown in Fig. S6 NPs respectively for different in Fig. 8. The percentage of stress by polymer chains rises with the strain at low = 1.5 while that by NPs declines. Meanwhile, that by polymer chains first increases and then decreases with the strain at high 2.0 while that by NPs has an opposite trend. In addition, the polymer chains bear a less percentage of stress with the while NPs bear a more one. This indicates that the effective stress transfer from polymer to NPs can happen at 2.0. Last, the number of voids is characterized with the strain for different in Fig. 9(a). It is found that the number of voids at strain ≈ 0.22 is significantly reduced with the increase of . The polymer-NP interface is reinforced with the increase of . Thus, as shown in Fig. 9(b), the probability of voids at the polymer-NP interface is reduced with the while that in the matrix enhanced. In total, the A beads on the surface of NPs can absorb the polymer chains due to the strong interaction, which act as the temporary cross-linked points. This can help to form a polymer-NP network, which can enhance both the maximum stress and the elongation. Thus, the fracture energy of PNCs shows a continuous increase with the . Currently, the fracture property of PNCs has been investigated, which can be tuned by many factors such as the shape, size and concentration of fillers, the polymer-filler interaction and so on. [3−7] The nucleation and the evolution process of nanovoids is a key factor in understanding the fracture properties of PNCs. Zhang et al. [12−15] first employed the SAXS to explore the formation of nanovoids in carbon black filled styrene-butadiene rubber. The scattering invariant is used to obtain the void volume fraction. The resolution size of nanovoids ranges from 20 nm to 30 nm by instruments. However, the size of most initial nanovoids is less than ten nanometers, which cannot be detected. Since a large increase in scattering intensity near the primary beam is rather an indication of voids than a detection, such estimations are phenomenological. Thus, the nanovoids are very difficult to be accurately recorded in experiments. In addition, our adopted coarsegrained model simulates a common polymer system whose nanoscale is much smaller than the experimental scale. The mapping rule from this model to the experiments can be referred to a recent work by Everaers and coworkers. [34] Meanwhile, NPs tend to form the aggregation because of their selfattractive interaction in experiments. Thus, the surface modification of NPs is an effective method for enhancing their compatibility with polymer chains. [28] Besides, the shear and tensile fields are employed to break the aggregation down. Unfortunately, once these fields are removed, NPs will gradually aggregate. As a result, the dispersion state of NPs does not reach the thermodynamic equilibrium state in experiments while they are in the thermodynamic equilibrium state in the simulation. In other works, the dispersion states of NPs between the experiment and simulation are not completely the same. The quantitative comparison between the experiments and our simulation should be careful. In summary, by analyzing the stress contribution of different components and the evolution process of initial nanovoids, this simulation work gives a clear understanding on how the functionalized NPs influence the fracture behavior of PNCs, which has not been reported.

CONCLUSIONS
In this work, we explored the cavitation and crazing process of polymer nanocomposites (PNCs) during the fracture process via a coarse-grained molecular dynamics simulation. Considering the fracture energy, the optimal fracture property is realized at the moderate NP functionalization degree while it continuously increases with the interaction between polymer and modified beads. The bond orientation degree is first analyzed with respect to the strain whose trend is consistent with the elongation. By quantifying the stress contribution, the stress contributed by polymer chains gradually declines with the or the while that by NPs rises. It is noted that most stress is borne by polymer chains at low or the . Interesting, the percentage of stress by polymer chains exhibits a first increase and then a decrease with the strain while that by NPs shows an opposite trend. Consistent with the stress-strain behavior, the number of voids first increases and then decreases with increasing the strain. This reflects the coalescence process of small voids into large ones. In particular, the nucleation of voids prefers to occur from the polymer-NP interface to the polymer matrix with increasing the or . As a result, the maximum number of voids first rises and then declines with the while it continuously decreases with the . In general, this work can provide a fundamental understanding of the cavitation and crazing process of PNCs during the fracture process by adopting the NP functionalization.

Electronic Supplementary Information
Electronic supplementary information (ESI) is available free of charge in the online version of this article at http://dx.doi.org/ 10.1007/s10118-020-2488-5.